Heterogeneous Treatment Effects of Intensive Glycemic Control on Kidney Microvascular Outcomes and Mortality in ACCORD
Vivek Charu, Jane W. Liang, Glenn M. Chertow, Zhuo Jun Li, Maria E. Montez-Rath, Pascal Geldsetzer, Ian H. de Boer, Lu Tian, and Manjula Kurella Tamura
2024-01-04
library(tidyverse)
library(zoo)
library(survminer)
library(survival)
library(survRM2)
library(foreach)
library(doParallel)
# Cluster object
cl = parallel::makeCluster(parallel::detectCores() - 2)
# Register parallel backend
registerDoParallel(cl)
# Set background to be white for all ggplots
theme_set(theme_classic())Data processing
We read in the data and extract the outcome and covariate variables of interest. Duplicate observations and observations where the treatment arm or variables needed to calculate the KFRE are missing are dropped.
# Read in datasets
# Raw data
raw_data = read.csv("/Users/jliang/Library/CloudStorage/Box-Box/Jane-Vivek/ACCORD HTE Analysis/Data/New data 26SEP2022/accordall.csv", header = TRUE)
# Additional outcomes
other_df = read.csv("/Users/jliang/Library/CloudStorage/Box-Box/Jane-Vivek/ACCORD HTE Analysis/Data/ACCORD_GLY_NEWOUTCOMES-2.csv", header = TRUE)
hypo_df = read.csv("/Users/jliang/Library/CloudStorage/Box-Box/RELATE\ CKD/Study\ Datasets/ACCORD/ACCORD_2017b_2\ 2/Main_Study/3-Data_Sets-Analysis/3a-Analysis_Data_Sets/csv/hypoglycemiatime1st.csv")
cvd_df = read.csv("/Users/jliang/Library/CloudStorage/Box-Box/RELATE\ CKD/Study\ Datasets/ACCORD/ACCORD_2017b_2\ 2/Main_Study/3-Data_Sets-Analysis/3a-Analysis_Data_Sets/csv/cvdoutcomes.csv")
# Concomitant medications
meds_df = read.csv("/Users/jliang/Library/CloudStorage/Box-Box/RELATE\ CKD/Study\ Datasets/ACCORD/ACCORD_2017b_2\ 2/Main_Study/3-Data_Sets-Analysis/3a-Analysis_Data_Sets/csv/concomitantmeds.csv")
# Classification groups for medication
meds_groups_df = read.csv("Copy of Meds_ACCORD.csv")[,1:2]
names(meds_groups_df) = c("med", "group")
# Subset for death outcomes
other_df = other_df %>%
select(maskid = id, ALLDEATH, alldeath_fu, CVDEATH, cvdeath_fu)
# Subset for 1st assisted hypoglycemic event
hypo_df = hypo_df %>%
mutate(hypoglycemia = 1 - censor_any,
hypoglycemia_fu = 365 * fuyrs_any) %>%
select(maskid = MaskID, hypoglycemia, hypoglycemia_fu)
# Subset for CVD primary outcome
cvd_df = cvd_df %>%
mutate(cvd_primary = 1 - censor_po,
cvd_primary_fu = 365 * fuyrs_po) %>%
select(maskid = MaskID, cvd_primary, cvd_primary_fu)
# Use friendlier variable names
meds_groups_df = meds_groups_df %>%
filter(group != "") %>%
mutate(group = recode(group, "anti-htn" = "anti_htn",
"chol lowering" = "chol_lowering",
"oral DM" = "oral_DM"))
# Convert data frame of classification groups to list
con_meds = sort(unique(meds_groups_df$group))
meds_groups_list = sapply(con_meds, function(x) {
meds_groups_df$med[meds_groups_df$group == x]
})
# Only consider baseline meds
meds_df = meds_df %>% filter(Visit == "BLR")
# Create a new variable for each con med classification group
for (med in con_meds) {
meds_df[[med]] =
rowSums(meds_df[,meds_groups_list[[med]]], na.rm = TRUE) > 0
}
# Merge in death outcomes
raw_data = merge(raw_data, other_df, by = "maskid", all.x = TRUE)
# Merge in 1st assisted hypoglycemic event
raw_data = merge(raw_data, hypo_df, by = "maskid", all.x = TRUE)
# Merge in CVD outcomes
raw_data = merge(raw_data, cvd_df, by = "maskid", all.x = TRUE)
# Merge in con meds
raw_data = merge(raw_data %>% select(!diuretic),
meds_df %>% select(maskid = MaskID, all_of(con_meds)),
by = "maskid", all.x = TRUE)
# Create composite kidney outcome
raw_data$neph_composite = pmax(raw_data$Neph2, raw_data$Neph3)
raw_data$neph_composite_fu = pmin(raw_data$Neph2Days, raw_data$Neph3Days)
# Quantitative variables to potentially include in the models
subset_quant = c("female", "baseline_age", "raceth",
"sbp", "dbp", "pulsepres",
"chol", "trig", "vldl", "ldl", "hdl",
"alt", "cpk",
"fpg", "hba1c",
"ualb", "ucreat", "uacr",
"ckd2021GFR", "screat",
"bmi",
"smokelif",
con_meds)
# List of outcome variables
outcome_list = list(
neph_composite = c(status = "neph_composite", time = "neph_composite_fu"),
Neph2 = c(status = "Neph2", time = "Neph2Days"),
Neph3 = c(status = "Neph3", time = "Neph3Days"),
ALLDEATH = c(status = "ALLDEATH", time = "alldeath_fu"),
CVDEATH = c(status = "CVDEATH", time = "cvdeath_fu"),
hypoglycemia = c(status = "hypoglycemia", time = "hypoglycemia_fu"),
cvd_primary = c(status = "cvd_primary", time = "cvd_primary_fu"))
# All outcomes
outcomes = unlist(outcome_list, use.names = FALSE)
# Race should be a factor
raw_data$raceth = as.factor(raw_data$raceth)
# Make relevant exclusions and define new dataset
df = raw_data %>%
select(all_of(c("maskid", "glyarm", "kfrs", outcomes, subset_quant))) %>%
drop_na(all_of(c("glyarm", "kfrs"))) %>% # Drop those with missing values
unique() # Drop duplicates
# Include race contrasts
df = data.frame(df,
raceth0 = ifelse(df$raceth == 0, 1, 0),
model.matrix(~raceth, df)[,-1])# Read in the longitudinal labs data and calculate eGFR.
dir = "/Users/jliang/Library/CloudStorage/Box-Box/RELATE\ CKD/Study\ Datasets/ACCORD/ACCORD_2017b_2\ 2/Main_Study/3-Data_Sets-Analysis/3a-Analysis_Data_Sets/csv/"
activitystatus = read.csv(paste0(dir, "activitystatus.csv"))
otherlabs = read.csv(paste0(dir, "otherlabs.csv"))
accord_key = read.csv(paste0(dir, "accord_key.csv"))
# Calculate eGFR
calc_eGFR = function(AGE, MALE, SCR) {
if (length(AGE) != length(MALE) ||
length(AGE) != length(SCR) ||
length(MALE) != length(SCR)) {
warning("Input arguments are not the same length.")
}
kau = 0.7*(1 - MALE) + 0.9*MALE
alpha = -0.241*(1 - MALE) - 0.302*MALE
beta = 1.012*(1 - MALE) + 1*MALE
return(142 * pmin(SCR/kau, 1)^alpha * pmax(SCR/kau, 1)^(-1.2) *
0.9938^AGE * beta)
}
# Combine otherlabs (contains serum creatinine) with activitystatus (contains
# FU time) and accord_key (contains age, sex, and arm)
# Calculate eGFR and 6-month spline
# Only include observations with 7 years and for IDs included in the main df
egfr_df = right_join(activitystatus %>% select(MaskID, Visit, days_from_baseline),
otherlabs %>% select(MaskID, Visit, screat, gfr),
by = c("MaskID", "Visit")) %>%
left_join(accord_key %>%
select(MaskID, female, baseline_age, raceclass, arm, treatment),
by = "MaskID") %>%
mutate(treat = ifelse(arm %in% c(3, 4, 7, 8), 1, 0),
baseline_age = baseline_age + days_from_baseline / 365) %>%
mutate(egfr = calc_eGFR(baseline_age, female == 0, screat),
spline = ifelse(days_from_baseline >= 180,
days_from_baseline - 180, 0)) %>%
filter(days_from_baseline <= 365*7,
MaskID %in% df$maskid) %>%
drop_na()
# Create data frame for time to reduction of eGFR by 40% of baseline
egfr40_df = left_join(
egfr_df %>%
group_by(MaskID) %>%
arrange(desc(days_from_baseline)) %>%
slice(1) %>%
ungroup() %>%
select(MaskID, days_from_baseline),
egfr_df %>%
group_by(MaskID) %>%
arrange(days_from_baseline) %>%
mutate(egfr40_cumsum = cumsum(egfr <= egfr[days_from_baseline == 0] * 0.6),
egfr40_consec_sum = rollsum(egfr <= egfr[days_from_baseline == 0] * 0.6,
k = 2, na.pad = TRUE, align = "right")) %>%
summarize(
# 40% eGFR reduction for 2 consecutive measures
egfr40_2consec = as.numeric(any(egfr40_consec_sum == 2, na.rm = TRUE)),
egfr40_2consec_fu = first(days_from_baseline[egfr40_consec_sum == 2],
na_rm = TRUE)) %>%
ungroup(),
by = "MaskID") %>%
rename(maskid = MaskID) %>%
mutate(egfr40_2consec = replace_na(egfr40_2consec, 0),
egfr40_2consec_fu = coalesce(egfr40_2consec_fu, days_from_baseline))
# Merge in event status and time to reduction of eGFR by 40% of baseline
df = left_join(df, egfr40_df, by = "maskid")
# Create composite out for first occurrence of 40% eGFR reduction or Neph3
# 40% eGFR reduction for 2 consecutive measures
df$egfr40_2consec = pmax(df$egfr40_2consec, df$Neph3)
df$egfr40_2consec_fu = pmin(df$egfr40_2consec_fu, df$Neph3Days)
# List of outcome variables
outcome_list = c(
outcome_list,
list(egfr40_2consec = c(status = "egfr40_2consec", time = "egfr40_2consec_fu")))
# All outcomes
outcomes = unlist(outcome_list, use.names = FALSE)We define a data frame that standardizes all continuous variables to have mean 0 and standard deviation 1.
# All covariates considered + race split into contrasts
expanded_covariates = c(setdiff(subset_quant, "raceth"),
paste0("raceth", 0:3))
# Binary variables
binary_vars = apply(df[,expanded_covariates], 2, function(x){
length(unique(x)) == 2 && all(sort(unique(x)) == c(0, 1))
})
# Data where continuous covariates are standardized
# Binary variables are left alone
std_df = sapply(1:length(expanded_covariates), function(i) {
if (binary_vars[i] == FALSE) {
return(scale(df[,expanded_covariates][,i]))
} else {
return(df[,expanded_covariates][,i])
}
})
colnames(std_df) = expanded_covariatesKidney failure risk equation (KFRE)
We calculate the 5-year predicted risk using the KFRE, define KFRE quartiles, and plot the distributions within treatment arms.
# Calculate 5-year KFRE
df$kfre5 = with(df, {
1 - 0.8996^exp(-0.2201 * (baseline_age/10 - 7.036) +
0.2467 * ((female==0) - 0.5642) -
0.5567 * (ckd2021GFR/5 - 7.222) +
0.4510 * ( log(uacr) - 5.137))
})
# kfre quartile thresholds
kfre_quart_thresh = quantile(df$kfre5, seq(0.25, 0.75, by = 0.25))
# Create groups based on thresholds
df$kfre_quarts = as.numeric(cut(df$kfre5,
breaks = c(-Inf, as.numeric(kfre_quart_thresh), Inf),
labels = c(1:(length(kfre_quart_thresh)+1))))
# Histograms of KFRE for treatment and control
df %>%
ggplot(aes(x = kfre5, fill = as.factor(glyarm), color = as.factor(glyarm))) +
geom_histogram(alpha = 0.3, position="identity") +
scale_x_log10(breaks = 10^(-6:0),
labels = function(x) paste0(formatC(x*100, format = "fg"), "%")) +
geom_vline(xintercept = kfre_quart_thresh, color = "gray60") +
geom_text(data = data.frame(kfre_quart_thresh = kfre_quart_thresh),
aes(x = kfre_quart_thresh, y = Inf,
label = sprintf("%1.3f%%", 100*kfre_quart_thresh)),
inherit.aes = FALSE, color = "grey60",
hjust = 1.1, vjust = 1.4, size = 3, angle = 90) +
xlab("5-year predicted risk by KFRE (log)") + ylab("") +
scale_fill_discrete(name = "", labels = c("0" = "Control", "1" = "Treatment")) +
scale_color_discrete(name = "", labels = c("0" = "Control", "1" = "Treatment"))
ggsave("figs_and_tabs/kfre_hist.png", width = 6, height = 4)
# General summary function that returns the mean, standard deviation, median
# and IQR of x
my_summary = function(x, na.rm = TRUE) {
c("Mean" = mean(x, na.rm = na.rm),
"SD" = sd(x, na.rm = na.rm),
"Median" = median(x, na.rm = na.rm),
"IQR" = IQR(x, na.rm = na.rm))
}
# Distribution of KFRE for treatment and control
kfre_dist_tab = data.frame(
do.call(rbind, c(list(my_summary(df$kfre5)),
tapply(df$kfre5, df$glyarm, my_summary))),
row.names = c("Overall", "Control", "Treatment")
)
write.csv(kfre_dist_tab %>% round(5), file = "figs_and_tabs/kfre_dist_tab.csv")
kfre_dist_tab %>% round(5)## Mean SD Median IQR
## Overall 0.00259 0.01875 0.00014 0.00075
## Control 0.00226 0.01229 0.00013 0.00071
## Treatment 0.00291 0.02346 0.00015 0.00079These are summary statistics for the scores within each KFRE quartile.
# Summary stats for KFRE within quartile
kfre_quart_tab = data.frame(do.call(rbind,
lapply(1:4, function(i) {
c(Quartile = i,
my_summary(df$kfre5[df$kfre_quarts==i]))
})
))
write.csv(kfre_quart_tab %>% round(5),
file = "figs_and_tabs/kfre_quart_tab.csv", row.names = FALSE)
kfre_quart_tab %>% round(5)## Quartile Mean SD Median IQR
## 1 1 0.00002 0.00001 0.00002 0.00001
## 2 2 0.00008 0.00003 0.00007 0.00005
## 3 3 0.00034 0.00017 0.00029 0.00026
## 4 4 0.00991 0.03653 0.00275 0.00575Number of observations overall and in each arm.
c(table(raw_data$glyarm), Overall = nrow(raw_data))## 0 1 Overall
## 5126 5140 10266Number of observations overall and in each arm, restricting to those with non-missing covariates needed to calculate KFRE.
c(table(df$glyarm), Overall = nrow(df))## 0 1 Overall
## 4873 4904 9777Number of non-missing observations for each outcome.
t(sapply(outcome_list, function(x) {
my_df = na.omit(df[,c(x["status"], x["time"], "glyarm")])
c(table(my_df$glyarm), Overall = nrow(my_df))
}))## 0 1 Overall
## neph_composite 4237 4219 8456
## Neph2 4237 4219 8456
## Neph3 4865 4895 9760
## ALLDEATH 4873 4904 9777
## CVDEATH 4873 4904 9777
## hypoglycemia 4873 4904 9777
## cvd_primary 4873 4904 9777
## egfr40_2consec 4865 4895 9760Descriptive statistics
For each quartile and overall, we report the proportion (standard deviation) for all binary variables and the median (IQR) of all continuous variables in the dataset.
# Data dictionary for nice variable names
dat_dict = data.frame(rbind(
c("female", "Female gender", "(yes/no)"),
c("raceth0", "White", "(yes/no)"),
c("raceth1", "Black", "(yes/no)"),
c("raceth3", "Hispanic", "(yes/no)"),
c("raceth2", "Other race/ethnicity", "(yes/no)"),
c("baseline_age", "Baseline age", "(years)"),
c("bmi", "Body mass index", "(kg/m2)"),
c("sbp", "Systolic blood pressure", "(mmHg)"),
c("dbp", "Diastolic blood pressure", "(mmHg)"),
c("pulsepres", "Pulse pressure", "(mmHg)"),
c("chol", "Total cholesterol", "(mg/dL)"),
c("trig", "Triglyercides", "(mg/dL)"),
c("vldl", "Very low density lipoprotein", "(mg/dL)"),
c("ldl", "Low density lipoprotein", "(mg/dL)"),
c("hdl", "High density lipoprotein", "(mg/dL)"),
c("alt", "Alanine aminotransferase", "(mg/dL)"),
c("cpk", "Creatine phosphokinase", "(mg/dL)"),
c("fpg", "Fasting plasma glucose", "(mg/dL)"),
c("hba1c", "Glycosylated hemoglobin", "(%)"),
c("ualb", "Urinary albumin", "(mg/dL)"),
c("ucreat", "Urinary creatinine", "(mg/dL)"),
c("uacr", "Urinary albumin to creatinine ratio", "(mg/g)"),
c("ckd2021GFR", "Estimated glomerular filtration rate", "(mL/min/1.73m2)"),
c("screat", "Serum creatinine", "(mg/dL)"),
c("anti_htn", "Anti-hypertension medication", "(yes/no)"),
c("chol_lowering", "Cholesterol-lowering medication", "(yes/no)"),
c("diuretic", "Diuretic medication", "(yes/no)"),
c("insulin", "Insulin", "(yes/no)"),
c("oral_DM", "Oral diabetes medication", "(yes/no)"),
c("neph_composite", "Composite kidney outcome", "(yes/no)"),
c("Neph2", "Development of macro-albuminuria", "(yes/no)"),
c("Neph3", "Renal failure or ESRD (dialysis) or SCr>3.3", "(yes/no)"),
c("ALLDEATH", "All cause death", "(yes/no)"),
c("CVDEATH", "Cardiovascular death", "(yes/no)"),
c("hypoglycemia", "1st assisted hypoglycemic event", "(yes/no)"),
c("egfr40_2consec", "ESRD or sustained (consecutive) eGFR reduction of 40%", "(yes/no)")
))
names(dat_dict) = c("short", "long", "units")
# Center (median or mean) of each variable
center_tab = cbind(
# Overall
sapply(c(expanded_covariates, names(outcome_list)), function(var) {
if (var %in% c("glyarm", names(outcome_list),
names(binary_vars)[binary_vars])) { # Binary variables
mean(df[,var], na.rm = TRUE)
} else { # Continuous
median(df[,var], na.rm = TRUE)
}
}),
# By quartile
sapply(1:4, function(i) {
sapply(c(expanded_covariates, names(outcome_list)), function(var) {
if (var %in% c("glyarm", names(outcome_list),
names(binary_vars)[binary_vars])) { # Binary variables
mean(df[df$kfre_quarts==i,var], na.rm = TRUE)
} else { # Continuous
median(df[df$kfre_quarts==i,var], na.rm = TRUE)
}
})
})
)
# Spread (IQR or standard deviation) of each variable
spread_tab = cbind(
# Overall
sapply(c(expanded_covariates, names(outcome_list)), function(var) {
if (var %in% c("glyarm", names(outcome_list),
names(binary_vars)[binary_vars])) { # Binary variables
sd(df[,var], na.rm = TRUE)
} else { # Continuous
IQR(df[,var], na.rm = TRUE)
}
}),
# By quartile
sapply(1:4, function(i){
sapply(c(expanded_covariates, names(outcome_list)), function(var) {
if (var %in% c("glyarm", names(outcome_list),
names(binary_vars)[binary_vars])) { # Binary variables
sd(df[df$kfre_quarts==i,var], na.rm = TRUE)
} else { # Continuous
IQR(df[df$kfre_quarts==i,var], na.rm = TRUE)
}
})
})
)
# Number of non-missing values for each covariate
N_tab = sapply(c(expanded_covariates, names(outcome_list)), function(var) {
sum(!is.na(df[,var]))
})
# Create Table 1
tab1 = matrix(paste0(round(center_tab, 2), " (", round(spread_tab, 2), ")"),
nrow = nrow(center_tab), ncol = 5)
rownames(tab1) = rownames(center_tab)
colnames(tab1) = c("Overall", paste("Quartile", 1:4))
tab1 = cbind(tab1, N = N_tab)
# Use more descriptive variable names
tab1 = merge(dat_dict, tab1, by.x = "short", by.y = 0, sort = FALSE)
rownames(tab1) = paste(tab1$long, tab1$units)
tab1$short = NULL; tab1$long = NULL; tab1$units = NULL
write.csv(tab1, file = "figs_and_tabs/tab1.csv")
tab1## Overall
## Female gender (yes/no) 0.38 (0.49)
## White (yes/no) 0.62 (0.48)
## Black (yes/no) 0.19 (0.39)
## Hispanic (yes/no) 0.07 (0.26)
## Other race/ethnicity (yes/no) 0.11 (0.32)
## Baseline age (years) 62.1 (9.5)
## Body mass index (kg/m2) 31.83 (7.71)
## Systolic blood pressure (mmHg) 136 (22)
## Diastolic blood pressure (mmHg) 75 (14)
## Pulse pressure (mmHg) 60 (19)
## Total cholesterol (mg/dL) 178 (53)
## Triglyercides (mg/dL) 156 (122)
## Very low density lipoprotein (mg/dL) 31 (25)
## Low density lipoprotein (mg/dL) 101 (44)
## High density lipoprotein (mg/dL) 40 (14)
## Alanine aminotransferase (mg/dL) 24 (14)
## Creatine phosphokinase (mg/dL) 106 (92)
## Fasting plasma glucose (mg/dL) 168 (65)
## Glycosylated hemoglobin (%) 8.1 (1.3)
## Urinary albumin (mg/dL) 1.69 (4.48)
## Urinary creatinine (mg/dL) 116.9 (77.5)
## Urinary albumin to creatinine ratio (mg/g) 14 (40)
## Estimated glomerular filtration rate (mL/min/1.73m2) 87.95 (26.47)
## Serum creatinine (mg/dL) 0.9 (0.2)
## Anti-hypertension medication (yes/no) 0.81 (0.39)
## Cholesterol-lowering medication (yes/no) 0.68 (0.47)
## Diuretic medication (yes/no) 0.37 (0.48)
## Insulin (yes/no) 0.35 (0.48)
## Oral diabetes medication (yes/no) 0.83 (0.37)
## Composite kidney outcome (yes/no) 0.1 (0.29)
## Development of macro-albuminuria (yes/no) 0.07 (0.26)
## Renal failure or ESRD (dialysis) or SCr>3.3 (yes/no) 0.03 (0.17)
## All cause death (yes/no) 0.07 (0.26)
## Cardiovascular death (yes/no) 0.03 (0.18)
## 1st assisted hypoglycemic event (yes/no) 0.12 (0.33)
## ESRD or sustained (consecutive) eGFR reduction of 40% (yes/no) 0.11 (0.31)
## Quartile 1
## Female gender (yes/no) 0.48 (0.5)
## White (yes/no) 0.66 (0.47)
## Black (yes/no) 0.11 (0.31)
## Hispanic (yes/no) 0.09 (0.28)
## Other race/ethnicity (yes/no) 0.14 (0.35)
## Baseline age (years) 58.6 (6)
## Body mass index (kg/m2) 32.27 (8.13)
## Systolic blood pressure (mmHg) 133 (20)
## Diastolic blood pressure (mmHg) 76 (13)
## Pulse pressure (mmHg) 56 (17)
## Total cholesterol (mg/dL) 182 (54)
## Triglyercides (mg/dL) 161 (131)
## Very low density lipoprotein (mg/dL) 32 (26)
## Low density lipoprotein (mg/dL) 102 (45)
## High density lipoprotein (mg/dL) 40 (13)
## Alanine aminotransferase (mg/dL) 25 (14)
## Creatine phosphokinase (mg/dL) 91 (74)
## Fasting plasma glucose (mg/dL) 173 (67)
## Glycosylated hemoglobin (%) 8.1 (1.3)
## Urinary albumin (mg/dL) 0.9 (1.09)
## Urinary creatinine (mg/dL) 112.7 (73.4)
## Urinary albumin to creatinine ratio (mg/g) 8 (9)
## Estimated glomerular filtration rate (mL/min/1.73m2) 102.37 (6.11)
## Serum creatinine (mg/dL) 0.7 (0.2)
## Anti-hypertension medication (yes/no) 0.74 (0.44)
## Cholesterol-lowering medication (yes/no) 0.65 (0.48)
## Diuretic medication (yes/no) 0.25 (0.44)
## Insulin (yes/no) 0.28 (0.45)
## Oral diabetes medication (yes/no) 0.85 (0.36)
## Composite kidney outcome (yes/no) 0.05 (0.21)
## Development of macro-albuminuria (yes/no) 0.02 (0.13)
## Renal failure or ESRD (dialysis) or SCr>3.3 (yes/no) 0.03 (0.16)
## All cause death (yes/no) 0.04 (0.19)
## Cardiovascular death (yes/no) 0.01 (0.12)
## 1st assisted hypoglycemic event (yes/no) 0.08 (0.27)
## ESRD or sustained (consecutive) eGFR reduction of 40% (yes/no) 0.08 (0.27)
## Quartile 2
## Female gender (yes/no) 0.3 (0.46)
## White (yes/no) 0.67 (0.47)
## Black (yes/no) 0.14 (0.35)
## Hispanic (yes/no) 0.07 (0.26)
## Other race/ethnicity (yes/no) 0.12 (0.33)
## Baseline age (years) 62 (8.6)
## Body mass index (kg/m2) 31.81 (7.47)
## Systolic blood pressure (mmHg) 135 (21)
## Diastolic blood pressure (mmHg) 75 (14)
## Pulse pressure (mmHg) 60 (18)
## Total cholesterol (mg/dL) 177 (53)
## Triglyercides (mg/dL) 156 (121.5)
## Very low density lipoprotein (mg/dL) 31 (25)
## Low density lipoprotein (mg/dL) 99 (43)
## High density lipoprotein (mg/dL) 39 (13)
## Alanine aminotransferase (mg/dL) 25 (15)
## Creatine phosphokinase (mg/dL) 102 (83)
## Fasting plasma glucose (mg/dL) 170 (63)
## Glycosylated hemoglobin (%) 8.1 (1.3)
## Urinary albumin (mg/dL) 1.94 (4.14)
## Urinary creatinine (mg/dL) 116.3 (80.1)
## Urinary albumin to creatinine ratio (mg/g) 17 (40)
## Estimated glomerular filtration rate (mL/min/1.73m2) 95.91 (9.16)
## Serum creatinine (mg/dL) 0.8 (0.2)
## Anti-hypertension medication (yes/no) 0.8 (0.4)
## Cholesterol-lowering medication (yes/no) 0.68 (0.47)
## Diuretic medication (yes/no) 0.32 (0.46)
## Insulin (yes/no) 0.3 (0.46)
## Oral diabetes medication (yes/no) 0.85 (0.35)
## Composite kidney outcome (yes/no) 0.09 (0.29)
## Development of macro-albuminuria (yes/no) 0.07 (0.25)
## Renal failure or ESRD (dialysis) or SCr>3.3 (yes/no) 0.02 (0.15)
## All cause death (yes/no) 0.06 (0.24)
## Cardiovascular death (yes/no) 0.02 (0.15)
## 1st assisted hypoglycemic event (yes/no) 0.09 (0.29)
## ESRD or sustained (consecutive) eGFR reduction of 40% (yes/no) 0.11 (0.31)
## Quartile 3
## Female gender (yes/no) 0.34 (0.47)
## White (yes/no) 0.61 (0.49)
## Black (yes/no) 0.23 (0.42)
## Hispanic (yes/no) 0.06 (0.25)
## Other race/ethnicity (yes/no) 0.1 (0.3)
## Baseline age (years) 63.6 (9.4)
## Body mass index (kg/m2) 31.7 (7.64)
## Systolic blood pressure (mmHg) 136 (23)
## Diastolic blood pressure (mmHg) 75 (15)
## Pulse pressure (mmHg) 61 (18)
## Total cholesterol (mg/dL) 176 (51)
## Triglyercides (mg/dL) 148 (115)
## Very low density lipoprotein (mg/dL) 30 (23)
## Low density lipoprotein (mg/dL) 100 (42)
## High density lipoprotein (mg/dL) 40 (13)
## Alanine aminotransferase (mg/dL) 24 (13.12)
## Creatine phosphokinase (mg/dL) 113 (102)
## Fasting plasma glucose (mg/dL) 167 (63.25)
## Glycosylated hemoglobin (%) 8.1 (1.2)
## Urinary albumin (mg/dL) 1.85 (5.54)
## Urinary creatinine (mg/dL) 121.3 (78.55)
## Urinary albumin to creatinine ratio (mg/g) 15 (44.5)
## Estimated glomerular filtration rate (mL/min/1.73m2) 81.83 (9.54)
## Serum creatinine (mg/dL) 1 (0.2)
## Anti-hypertension medication (yes/no) 0.83 (0.37)
## Cholesterol-lowering medication (yes/no) 0.7 (0.46)
## Diuretic medication (yes/no) 0.39 (0.49)
## Insulin (yes/no) 0.37 (0.48)
## Oral diabetes medication (yes/no) 0.84 (0.37)
## Composite kidney outcome (yes/no) 0.09 (0.29)
## Development of macro-albuminuria (yes/no) 0.07 (0.26)
## Renal failure or ESRD (dialysis) or SCr>3.3 (yes/no) 0.02 (0.15)
## All cause death (yes/no) 0.07 (0.25)
## Cardiovascular death (yes/no) 0.03 (0.18)
## 1st assisted hypoglycemic event (yes/no) 0.13 (0.34)
## ESRD or sustained (consecutive) eGFR reduction of 40% (yes/no) 0.11 (0.32)
## Quartile 4
## Female gender (yes/no) 0.39 (0.49)
## White (yes/no) 0.57 (0.5)
## Black (yes/no) 0.29 (0.45)
## Hispanic (yes/no) 0.06 (0.24)
## Other race/ethnicity (yes/no) 0.09 (0.28)
## Baseline age (years) 65.15 (10.3)
## Body mass index (kg/m2) 31.54 (7.6)
## Systolic blood pressure (mmHg) 138 (23)
## Diastolic blood pressure (mmHg) 73 (16)
## Pulse pressure (mmHg) 65 (21)
## Total cholesterol (mg/dL) 178 (52)
## Triglyercides (mg/dL) 158 (121)
## Very low density lipoprotein (mg/dL) 32 (24)
## Low density lipoprotein (mg/dL) 101 (44.5)
## High density lipoprotein (mg/dL) 40 (15)
## Alanine aminotransferase (mg/dL) 22 (12)
## Creatine phosphokinase (mg/dL) 120 (113)
## Fasting plasma glucose (mg/dL) 163 (67)
## Glycosylated hemoglobin (%) 8.2 (1.3)
## Urinary albumin (mg/dL) 3.99 (15.82)
## Urinary creatinine (mg/dL) 116.75 (76.23)
## Urinary albumin to creatinine ratio (mg/g) 35 (153)
## Estimated glomerular filtration rate (mL/min/1.73m2) 63.3 (13.24)
## Serum creatinine (mg/dL) 1.2 (0.3)
## Anti-hypertension medication (yes/no) 0.87 (0.33)
## Cholesterol-lowering medication (yes/no) 0.71 (0.46)
## Diuretic medication (yes/no) 0.51 (0.5)
## Insulin (yes/no) 0.45 (0.5)
## Oral diabetes medication (yes/no) 0.78 (0.41)
## Composite kidney outcome (yes/no) 0.17 (0.37)
## Development of macro-albuminuria (yes/no) 0.15 (0.36)
## Renal failure or ESRD (dialysis) or SCr>3.3 (yes/no) 0.04 (0.2)
## All cause death (yes/no) 0.12 (0.32)
## Cardiovascular death (yes/no) 0.06 (0.24)
## 1st assisted hypoglycemic event (yes/no) 0.18 (0.39)
## ESRD or sustained (consecutive) eGFR reduction of 40% (yes/no) 0.12 (0.33)
## N
## Female gender (yes/no) 9777
## White (yes/no) 9777
## Black (yes/no) 9777
## Hispanic (yes/no) 9777
## Other race/ethnicity (yes/no) 9777
## Baseline age (years) 9777
## Body mass index (kg/m2) 9771
## Systolic blood pressure (mmHg) 9733
## Diastolic blood pressure (mmHg) 9733
## Pulse pressure (mmHg) 9733
## Total cholesterol (mg/dL) 9775
## Triglyercides (mg/dL) 9775
## Very low density lipoprotein (mg/dL) 9775
## Low density lipoprotein (mg/dL) 9775
## High density lipoprotein (mg/dL) 9776
## Alanine aminotransferase (mg/dL) 9777
## Creatine phosphokinase (mg/dL) 9776
## Fasting plasma glucose (mg/dL) 9777
## Glycosylated hemoglobin (%) 9763
## Urinary albumin (mg/dL) 9777
## Urinary creatinine (mg/dL) 9777
## Urinary albumin to creatinine ratio (mg/g) 9777
## Estimated glomerular filtration rate (mL/min/1.73m2) 9777
## Serum creatinine (mg/dL) 9777
## Anti-hypertension medication (yes/no) 9777
## Cholesterol-lowering medication (yes/no) 9777
## Diuretic medication (yes/no) 9777
## Insulin (yes/no) 9777
## Oral diabetes medication (yes/no) 9777
## Composite kidney outcome (yes/no) 8456
## Development of macro-albuminuria (yes/no) 8456
## Renal failure or ESRD (dialysis) or SCr>3.3 (yes/no) 9777
## All cause death (yes/no) 9777
## Cardiovascular death (yes/no) 9777
## 1st assisted hypoglycemic event (yes/no) 9777
## ESRD or sustained (consecutive) eGFR reduction of 40% (yes/no) 9777Balance between treatment and control arm within quartile
To check the balance between the treatment and control groups, we take the difference in standardized means for each covariate, within each treatment arm and quartile.
# Mean of each covariate across individuals within each quartile and treatment
# group, for kfirs
# Standardized means were used for continuous covariates
std_means_balance_kfre = lapply(expanded_covariates, function(var) {
tapply(std_df[,var], list(df$kfre_quarts, df$glyarm),
mean, na.rm = TRUE)
})
names(std_means_balance_kfre) = expanded_covariates
# Split into treatment and control
# Treatment
std_means_trt_kfre = sapply(expanded_covariates, function(var){
std_means_balance_kfre[[var]][,"1"]
}) %>% t()
# Control
std_means_control_kfre = sapply(expanded_covariates, function(var){
std_means_balance_kfre[[var]][,"0"]
}) %>% t()
# Difference between treatment and control
std_means_diff_kfre =
merge(dat_dict,
sapply(expanded_covariates, function(var){
std_means_balance_kfre[[var]][,"1"] -
std_means_balance_kfre[[var]][,"0"]
}) %>% t(),
by.x = "short", by.y = 0, sort = FALSE)
rownames(std_means_diff_kfre) = std_means_diff_kfre$long
std_means_diff_kfre = std_means_diff_kfre[,-c(1:3)]We plot the mean difference between treatment and control for each variable, within each quartile.
data.frame(est = gather(data.frame(std_means_diff_kfre))$value,
Quartile = rep(as.factor(1:4), each = nrow(std_means_diff_kfre))) %>%
mutate(y = rep(1:nrow(std_means_diff_kfre), 4) +
rep(seq(-0.2, 0.2, length = 4), each = nrow(std_means_diff_kfre))) %>%
ggplot(aes(x = as.numeric(est), y = y, color = Quartile)) +
geom_point() +
scale_color_manual(values = c("navy", "dodgerblue", "mediumpurple","firebrick")) +
scale_y_continuous(breaks = 1:nrow(std_means_diff_kfre),
labels = rownames(std_means_diff_kfre),
expand = c(0.02, 0.02)) +
xlab("Difference in standardized mean between treatment and control") + ylab("") +
ggtitle("Imbalance by KFRE quartile")
ggsave("figs_and_tabs/kfre_imbalance.png", width = 6, height = 5)
write.csv(std_means_diff_kfre %>% round(3),
file = "figs_and_tabs/std_means_diff_kfre.csv")
std_means_diff_kfre %>% round(3)## 1 2 3 4
## Female gender -0.019 0.006 0.027 0.010
## White -0.010 0.010 0.011 -0.001
## Black 0.016 0.002 -0.006 0.004
## Hispanic 0.009 -0.013 -0.013 -0.002
## Other race/ethnicity -0.016 0.002 0.008 -0.001
## Baseline age -0.022 0.036 -0.059 0.005
## Body mass index 0.003 -0.057 0.021 -0.006
## Systolic blood pressure -0.049 0.034 -0.012 -0.028
## Diastolic blood pressure -0.025 -0.035 0.005 -0.023
## Pulse pressure -0.040 0.066 -0.017 -0.016
## Total cholesterol 0.066 -0.041 0.005 -0.004
## Triglyercides 0.021 -0.011 0.051 0.005
## Very low density lipoprotein 0.030 -0.009 0.050 -0.011
## Low density lipoprotein 0.064 -0.044 -0.020 -0.001
## High density lipoprotein -0.014 0.000 -0.027 0.011
## Alanine aminotransferase 0.021 -0.036 0.009 -0.012
## Creatine phosphokinase 0.010 -0.032 -0.081 0.048
## Fasting plasma glucose 0.083 0.004 -0.064 -0.086
## Glycosylated hemoglobin 0.089 -0.024 -0.076 -0.056
## Urinary albumin -0.001 -0.003 -0.022 -0.021
## Urinary creatinine -0.064 -0.013 -0.014 0.098
## Urinary albumin to creatinine ratio 0.005 -0.005 -0.006 0.037
## Estimated glomerular filtration rate 0.007 0.003 0.015 -0.026
## Serum creatinine 0.008 -0.015 -0.025 0.039
## Anti-hypertension medication -0.031 -0.018 0.006 -0.014
## Cholesterol-lowering medication -0.027 0.002 0.015 -0.006
## Diuretic medication -0.010 -0.006 -0.028 0.021
## Insulin -0.017 -0.006 -0.029 -0.021
## Oral diabetes medication -0.009 -0.009 0.019 0.010No covariates were imbalanced (defined as the mean difference between treatment and control being greater than 0.1 for at least one quartile), so we will not run an adjusted analysis.
# Which variables are imbalanced by more than 0.05 standard units?
is_imbalanced = apply(std_means_diff_kfre, 1, function(x){
any(abs(x) > 0.1)
})
imbalanced_vars = dat_dict$short[which(dat_dict$long %in%
names(is_imbalanced)[is_imbalanced])]
imbalanced_vars## character(0)Treatment effects
We estimate the 7-year RMST differences for each subgroup defined by the KFRE quartiles.
# Horizon
horizon = 365*7
# RMST difference overall
overall_rmst_fits = lapply(outcome_list, function(x){
df_sub = na.omit(df[,c(x["time"], x["status"], "glyarm", "kfre_quarts")])
if (x["status"] == "egfr40_2consec") {
horizon = 2525
} else {
horizon = 365*7
}
rmst2(df_sub[,x["time"]],
df_sub[,x["status"]],
df_sub[,"glyarm"],
tau = horizon)
})
names(overall_rmst_fits) = names(outcome_list)
overall_rmst = sapply(names(outcome_list), function(outcome){
overall_rmst_fits[[outcome]]$unadjusted.result[1,1]
})
# RMST difference by quartile
kfre_quart_rmst_fits = lapply(outcome_list, function(x){
df_sub = na.omit(df[,c(x["time"], x["status"], "glyarm", "kfre_quarts")])
lapply(1:4, function(i){
rmst2(df_sub[df_sub$kfre_quarts==i,x["time"]],
df_sub[df_sub$kfre_quarts==i,x["status"]],
df_sub[df_sub$kfre_quarts==i,"glyarm"],
tau = horizon)
})
})
names(kfre_quart_rmst_fits) = names(outcome_list)
kfre_quart_rmst = sapply(names(outcome_list), function(outcome){
sapply(1:4, function(i){
kfre_quart_rmst_fits[[outcome]][[i]]$unadjusted.result[1,1]
})
})We use 1000 bootstrap samples to obtain confidence intervals for the RMST differences. Tables of the RMST differences and normalized RMST differences (where the overall RMST difference is subtracted from each quartile’s estimate) with 95% bootstrap CIs are shown below.
# Number of bootstraps/permutations
B = 1000boot_kfre_rmst = foreach::foreach(
b = 1:B,
.packages = c("survRM2")
) %dopar% {
# Set seed
set.seed(b)
# Create bootstrap sample
boot_idx = sample(nrow(df), replace= TRUE)
df_boot = df[boot_idx,]
# Calculate overall RMST difference
overall_rmst = sapply(outcome_list, function(x){
df_sub = na.omit(df_boot[,c(x["time"], x["status"], "glyarm", "kfre_quarts")])
rmst2(df_sub[,x["time"]],
df_sub[,x["status"]],
df_sub[,"glyarm"],
tau = horizon)$unadjusted.result[1,1]
})
# Calculate RMST difference for each quartile
kfre_quart_rmst = sapply(outcome_list, function(x){
df_sub = na.omit(df_boot[,c(x["time"], x["status"], "glyarm", "kfre_quarts")])
if (x["status"] == "egfr40_2consec") {
horizon = 2525
} else {
horizon = 365*7
}
sapply(1:4, function(i){
rmst2(df_sub[df_sub$kfre_quarts==i,x["time"]],
df_sub[df_sub$kfre_quarts==i,x["status"]],
df_sub[df_sub$kfre_quarts==i,"glyarm"],
tau = horizon)$unadjusted.result[1,1]
})
})
return(list(overall_rmst = overall_rmst,
kfre_quart_rmst = kfre_quart_rmst))
}
save(boot_kfre_rmst, file = "boot_kfre_rmst.rData")# Overall
kfre_overall_rmst_boot_ci = apply(
sapply(boot_kfre_rmst, function(x){
x$overall_rmst
}), 1, quantile, c(0.025, 0.975))
# Not normalized
kfre_quart_rmst_boot_ci = lapply(1:4, function(i){
dat = sapply(boot_kfre_rmst, function(x){
x$kfre_quart_rmst[i,]
})
apply(dat, 1, quantile, c(0.025, 0.975))
})
# Make tables
kfre_quart_rmst_boot_tabs = lapply(names(outcome_list), function(outcome){
out = rbind(
data.frame(Quartile = "Overall",
"Est" = overall_rmst[outcome],
t(kfre_overall_rmst_boot_ci[,outcome])),
data.frame(Quartile = 1:4,
"Est" = kfre_quart_rmst[,outcome],
t(sapply(kfre_quart_rmst_boot_ci, function(x){x[,outcome]}))))
rownames(out) = NULL
return(out)
})
names(kfre_quart_rmst_boot_tabs) = names(outcome_list)
# Format table for printing
rmst_tab = data.frame(
Outcome = rep(names(kfre_quart_rmst_boot_tabs),
each = nrow(kfre_quart_rmst_boot_tabs[[1]])),
do.call(rbind, kfre_quart_rmst_boot_tabs),
row.names = NULL
)
rmst_tab = merge(dat_dict[,c("short", "long")], rmst_tab,
by.x = "short", by.y = "Outcome", all.y = TRUE,
sort = FALSE)[,-1]
names(rmst_tab) = c("Outcome", "Quartile", "Est", "2.5%", "97.5%")
rmst_tab[,c("Est", "2.5%", "97.5%")] =
round(rmst_tab[,c("Est", "2.5%", "97.5%")], 2)
write.csv(rmst_tab,
file = "figs_and_tabs/rmst_tab.csv", row.names = FALSE)
rmst_tab## Outcome Quartile Est
## 1 Composite kidney outcome 1 10.14
## 2 Composite kidney outcome 2 40.07
## 3 Composite kidney outcome 3 50.08
## 4 Composite kidney outcome Overall 48.40
## 5 Composite kidney outcome 4 114.77
## 6 Development of macro-albuminuria Overall 42.49
## 7 Development of macro-albuminuria 1 10.99
## 8 Development of macro-albuminuria 2 31.32
## 9 Development of macro-albuminuria 4 116.49
## 10 Development of macro-albuminuria 3 33.46
## 11 Renal failure or ESRD (dialysis) or SCr>3.3 Overall 5.43
## 12 Renal failure or ESRD (dialysis) or SCr>3.3 1 -4.10
## 13 Renal failure or ESRD (dialysis) or SCr>3.3 2 7.53
## 14 Renal failure or ESRD (dialysis) or SCr>3.3 4 0.99
## 15 Renal failure or ESRD (dialysis) or SCr>3.3 3 18.22
## 16 All cause death Overall -23.64
## 17 All cause death 1 -21.03
## 18 All cause death 2 -19.53
## 19 All cause death 3 11.06
## 20 All cause death 4 -56.70
## 21 Cardiovascular death Overall -15.43
## 22 Cardiovascular death 2 -16.46
## 23 Cardiovascular death 3 11.88
## 24 Cardiovascular death 4 -35.01
## 25 Cardiovascular death 1 -17.82
## 26 1st assisted hypoglycemic event Overall -232.97
## 27 1st assisted hypoglycemic event 2 -216.06
## 28 1st assisted hypoglycemic event 3 -247.08
## 29 1st assisted hypoglycemic event 4 -288.38
## 30 1st assisted hypoglycemic event 1 -171.59
## 31 ESRD or sustained (consecutive) eGFR reduction of 40% Overall -1.86
## 32 ESRD or sustained (consecutive) eGFR reduction of 40% 1 19.97
## 33 ESRD or sustained (consecutive) eGFR reduction of 40% 2 -4.53
## 34 ESRD or sustained (consecutive) eGFR reduction of 40% 3 12.84
## 35 ESRD or sustained (consecutive) eGFR reduction of 40% 4 -29.32
## 36 <NA> Overall 13.46
## 37 <NA> 1 -9.14
## 38 <NA> 2 6.54
## 39 <NA> 3 41.38
## 40 <NA> 4 23.27
## 2.5% 97.5%
## 1 -23.90 43.13
## 2 -4.65 84.70
## 3 7.95 92.69
## 4 25.34 69.60
## 5 58.14 176.41
## 6 22.52 61.31
## 7 -8.12 30.61
## 8 -4.86 68.63
## 9 62.97 173.51
## 10 -6.22 71.98
## 11 -7.02 16.95
## 12 -30.36 19.86
## 13 -14.43 31.42
## 14 -27.60 27.66
## 15 -3.14 39.30
## 16 -42.20 -6.64
## 17 -47.70 6.66
## 18 -53.51 12.56
## 19 -20.23 42.14
## 20 -100.19 -17.50
## 21 -28.26 -3.57
## 22 -40.91 4.34
## 23 -12.68 39.15
## 24 -65.91 -5.41
## 25 -37.08 0.40
## 26 -259.47 -208.03
## 27 -263.96 -168.89
## 28 -302.38 -192.26
## 29 -346.09 -227.82
## 30 -217.89 -127.53
## 31 -25.85 19.87
## 32 -25.65 66.44
## 33 -52.13 41.24
## 34 -30.46 56.44
## 35 -75.96 16.64
## 36 -10.10 34.64
## 37 -42.10 26.45
## 38 -37.73 48.36
## 39 0.72 86.08
## 40 -33.49 77.24# Normalized (subtract overall RMST difference)
kfre_quart_rmst_norm_boot_ci = lapply(1:4, function(i){
dat = sapply(boot_kfre_rmst, function(x){
x$kfre_quart_rmst[i,] - x$overall_rmst
})
apply(dat, 1, quantile, c(0.025, 0.975))
})
# Make tables
kfre_quart_rmst_boot_norm_tabs = lapply(names(outcome_list), function(outcome){
data.frame(Quartile = 1:4,
"Est" = kfre_quart_rmst[,outcome] - overall_rmst[outcome],
t(sapply(kfre_quart_rmst_norm_boot_ci, function(x){x[,outcome]})))
})
names(kfre_quart_rmst_boot_norm_tabs) = names(outcome_list)
# Format table for printing
rmst_norm_tab = data.frame(
Outcome = rep(names(kfre_quart_rmst_boot_norm_tabs),
each = nrow(kfre_quart_rmst_boot_norm_tabs[[1]])),
do.call(rbind, kfre_quart_rmst_boot_norm_tabs),
row.names = NULL
)
rmst_norm_tab = merge(dat_dict[,c("short", "long")], rmst_norm_tab,
by.x = "short", by.y = "Outcome", all.y = TRUE,
sort = FALSE)[,-1]
names(rmst_norm_tab) = c("Outcome", "Quartile", "Est", "2.5%", "97.5%")
rmst_norm_tab[,c("Est", "2.5%", "97.5%")] =
round(rmst_norm_tab[,c("Est", "2.5%", "97.5%")], 2)
write.csv(rmst_norm_tab,
file = "figs_and_tabs/rmst_norm_tab.csv", row.names = FALSE)
rmst_norm_tab## Outcome Quartile Est
## 1 Composite kidney outcome 1 -38.25
## 2 Composite kidney outcome 2 -8.33
## 3 Composite kidney outcome 3 1.68
## 4 Composite kidney outcome 4 66.37
## 5 Development of macro-albuminuria 1 -31.50
## 6 Development of macro-albuminuria 2 -11.17
## 7 Development of macro-albuminuria 3 -9.02
## 8 Development of macro-albuminuria 4 74.00
## 9 Renal failure or ESRD (dialysis) or SCr>3.3 1 -9.53
## 10 Renal failure or ESRD (dialysis) or SCr>3.3 2 2.09
## 11 Renal failure or ESRD (dialysis) or SCr>3.3 3 12.79
## 12 Renal failure or ESRD (dialysis) or SCr>3.3 4 -4.45
## 13 All cause death 1 2.61
## 14 All cause death 2 4.11
## 15 All cause death 3 34.70
## 16 All cause death 4 -33.06
## 17 Cardiovascular death 1 -2.39
## 18 Cardiovascular death 2 -1.03
## 19 Cardiovascular death 3 27.31
## 20 Cardiovascular death 4 -19.58
## 21 1st assisted hypoglycemic event 1 61.38
## 22 1st assisted hypoglycemic event 2 16.92
## 23 1st assisted hypoglycemic event 3 -14.11
## 24 1st assisted hypoglycemic event 4 -55.41
## 25 ESRD or sustained (consecutive) eGFR reduction of 40% 1 21.82
## 26 ESRD or sustained (consecutive) eGFR reduction of 40% 2 -2.67
## 27 ESRD or sustained (consecutive) eGFR reduction of 40% 3 14.70
## 28 ESRD or sustained (consecutive) eGFR reduction of 40% 4 -27.46
## 29 <NA> 1 -22.61
## 30 <NA> 2 -6.93
## 31 <NA> 3 27.92
## 32 <NA> 4 9.81
## 2.5% 97.5%
## 1 -70.97 -6.26
## 2 -46.88 31.06
## 3 -35.00 36.96
## 4 19.51 118.09
## 5 -55.25 -7.48
## 6 -46.01 21.76
## 7 -42.38 23.75
## 8 28.96 122.59
## 9 -32.56 13.00
## 10 -17.36 24.02
## 11 -7.34 33.15
## 12 -27.63 18.20
## 13 -21.95 30.27
## 14 -27.23 32.29
## 15 7.29 63.81
## 16 -66.43 -0.53
## 17 -20.93 17.13
## 18 -23.66 19.05
## 19 5.59 50.96
## 20 -44.91 4.62
## 21 18.55 103.79
## 22 -28.47 57.71
## 23 -60.68 33.34
## 24 -104.52 -5.58
## 25 -21.43 60.40
## 26 -41.31 38.06
## 27 -24.19 55.58
## 28 -68.53 12.52
## 29 -53.37 9.51
## 30 -44.85 30.90
## 31 -8.48 64.56
## 32 -34.86 54.22These are plots of the RMST differences by outcome, with a horizontal reference line drawn at the overall/ATE point estimate. Shading denotes the ATE CI.
for (outcome in names(kfre_quart_rmst_boot_tabs)) {
ate_tab = kfre_quart_rmst_boot_tabs[[outcome]][1,]
p1 = kfre_quart_rmst_boot_tabs[[outcome]][-1,] %>%
mutate(x = 1:4) %>%
ggplot(aes(x = x, y = as.numeric(Est))) +
geom_pointrange(aes(ymin = X2.5., ymax = X97.5.), size = 0.25) +
geom_hline(yintercept = ate_tab[,"Est"], color = "grey") +
annotate("rect", xmin = -Inf, xmax = Inf,
ymin = ate_tab[,"X2.5."], ymax = ate_tab[,"X97.5."],
alpha = 0.2) +
scale_x_continuous(breaks = 1:4, labels = 1:4) +
theme(legend.position = 'bottom', legend.box = 'vertical') +
xlab("Quartile") + ylab("RMST difference") +
ggtitle(dat_dict$long[which(dat_dict$short == outcome)])
print(p1)
ggsave(paste0("figs_and_tabs/rmst_", outcome, ".png"),
width = 6, height = 4)
}
Normalized RMST differences.
for (outcome in names(kfre_quart_rmst_boot_norm_tabs)) {
p1 = kfre_quart_rmst_boot_norm_tabs[[outcome]] %>%
mutate(x = 1:4) %>%
ggplot(aes(x = x, y = as.numeric(Est))) +
geom_pointrange(aes(ymin = X2.5., ymax = X97.5.), size = 0.25) +
scale_x_continuous(breaks = 1:4, labels = 1:4) +
theme(legend.position = 'bottom', legend.box = 'vertical') +
xlab("Quartile") + ylab("Normalized RMST difference") +
geom_hline(yintercept = 0, color = "grey") +
ggtitle(dat_dict$long[which(dat_dict$short == outcome)])
print(p1)
ggsave(paste0("figs_and_tabs/rmst_norm_", outcome, ".png"),
width = 6, height = 4)
}
Cox PH models and cumulative incidence curves
For each outcome, we fit a Cox PH model to treatment for the entire dataset and each group defined by the KFRE quartiles. Estimates of the treatment hazard ratios are shown below.
# Fit Cox PH model for each outcome to glyarm
coxph_fits = lapply(names(outcome_list), function(outcome){
# Extract data of interest
x = outcome_list[[outcome]]
df_sub = na.omit(df[,c(x["time"], x["status"], "glyarm")])
# Formula for survfit object
form = formula(paste0("Surv(", x["time"], ", " ,
x["status"], ") ~ glyarm"))
# Fit survfit object
fit = coxph(form, data = df_sub)
# Manually insert formula
fit$call$formula = form
return(fit)
})
names(coxph_fits) = names(outcome_list)
# Fit Cox PH models within quartile
coxph_fits_by_quart = lapply(names(outcome_list), function(outcome){
# Extract data of interest
x = outcome_list[[outcome]]
df_sub = na.omit(df[,c(x["time"], x["status"], "glyarm", "kfre_quarts")])
# Formula for survfit object
form = formula(paste0("Surv(", x["time"], ", " ,
x["status"], ") ~ glyarm"))
all_fits = lapply(1:4, function(i){
# Subset data to quartile
df_sub = df_sub[df_sub$kfre_quarts==i, ]
# Fit survfit object
fit = coxph(form, data = df_sub)
# Manually insert formula
fit$call$formula = form
return(fit)
})
return(all_fits)
})
names(coxph_fits_by_quart) = names(outcome_list)
# Make tables
coxph_coeff_tabs = lapply(names(outcome_list), function(outcome){
out = rbind(
data.frame(Quartile = "Overall",
t(exp(c(coef(coxph_fits[[outcome]]), confint(coxph_fits[[outcome]]))))
),
data.frame(Quartile = 1:4,
do.call(rbind,
lapply(coxph_fits_by_quart[[outcome]], function(x){
exp(c(coef(x), confint(x)))
}))
))
rownames(out) = NULL
colnames(out) = c("Quartile", "Est", "X2.5.", "X97.5.")
return(out)
})
names(coxph_coeff_tabs) = names(outcome_list)
# Format table for printing
coeff_tab = data.frame(
Outcome = rep(names(coxph_coeff_tabs),
each = nrow(coxph_coeff_tabs[[1]])),
do.call(rbind, coxph_coeff_tabs),
row.names = NULL
)
coeff_tab = merge(dat_dict[,c("short", "long")], coeff_tab,
by.x = "short", by.y = "Outcome", all.y = TRUE,
sort = FALSE)[,-1]
names(coeff_tab) = c("Outcome", "Quartile", "Est", "2.5%", "97.5%")
coeff_tab[,c("Est", "2.5%", "97.5%")] =
round(coeff_tab[,c("Est", "2.5%", "97.5%")], 2)
write.csv(coeff_tab,
file = "figs_and_tabs/coeff_tab.csv", row.names = FALSE)
coeff_tab## Outcome Quartile Est 2.5%
## 1 Composite kidney outcome 1 0.93 0.63
## 2 Composite kidney outcome 2 0.78 0.59
## 3 Composite kidney outcome 3 0.73 0.54
## 4 Composite kidney outcome Overall 0.75 0.65
## 5 Composite kidney outcome 4 0.66 0.53
## 6 Development of macro-albuminuria Overall 0.72 0.61
## 7 Development of macro-albuminuria 1 0.73 0.39
## 8 Development of macro-albuminuria 2 0.77 0.56
## 9 Development of macro-albuminuria 4 0.62 0.49
## 10 Development of macro-albuminuria 3 0.77 0.56
## 11 Renal failure or ESRD (dialysis) or SCr>3.3 Overall 0.93 0.73
## 12 Renal failure or ESRD (dialysis) or SCr>3.3 1 1.09 0.68
## 13 Renal failure or ESRD (dialysis) or SCr>3.3 2 0.89 0.52
## 14 Renal failure or ESRD (dialysis) or SCr>3.3 4 1.03 0.69
## 15 Renal failure or ESRD (dialysis) or SCr>3.3 3 0.65 0.38
## 16 All cause death Overall 1.20 1.04
## 17 All cause death 1 1.37 0.91
## 18 All cause death 2 1.15 0.84
## 19 All cause death 3 0.95 0.70
## 20 All cause death 4 1.30 1.03
## 21 Cardiovascular death Overall 1.30 1.05
## 22 Cardiovascular death 2 1.48 0.89
## 23 Cardiovascular death 3 0.85 0.55
## 24 Cardiovascular death 4 1.40 1.01
## 25 Cardiovascular death 1 1.75 0.88
## 26 1st assisted hypoglycemic event Overall 2.91 2.56
## 27 1st assisted hypoglycemic event 2 3.90 2.85
## 28 1st assisted hypoglycemic event 3 2.99 2.32
## 29 1st assisted hypoglycemic event 4 2.41 1.97
## 30 1st assisted hypoglycemic event 1 3.05 2.22
## 31 ESRD or sustained (consecutive) eGFR reduction of 40% Overall 1.02 0.90
## 32 ESRD or sustained (consecutive) eGFR reduction of 40% 1 0.92 0.70
## 33 ESRD or sustained (consecutive) eGFR reduction of 40% 2 1.03 0.81
## 34 ESRD or sustained (consecutive) eGFR reduction of 40% 3 0.92 0.73
## 35 ESRD or sustained (consecutive) eGFR reduction of 40% 4 1.16 0.93
## 36 <NA> Overall 0.91 0.81
## 37 <NA> 1 1.08 0.77
## 38 <NA> 2 0.94 0.73
## 39 <NA> 3 0.80 0.62
## 40 <NA> 4 0.89 0.73
## 97.5%
## 1 1.37
## 2 1.04
## 3 0.97
## 4 0.86
## 5 0.83
## 6 0.84
## 7 1.39
## 8 1.06
## 9 0.79
## 10 1.07
## 11 1.18
## 12 1.76
## 13 1.51
## 14 1.53
## 15 1.10
## 16 1.40
## 17 2.07
## 18 1.58
## 19 1.29
## 20 1.64
## 21 1.62
## 22 2.48
## 23 1.31
## 24 1.95
## 25 3.47
## 26 3.30
## 27 5.35
## 28 3.84
## 29 2.94
## 30 4.17
## 31 1.15
## 32 1.22
## 33 1.31
## 34 1.17
## 35 1.46
## 36 1.03
## 37 1.51
## 38 1.23
## 39 1.01
## 40 1.09We plot cumulative incidence curves and risk tables for each outcome.
# Fit Kaplan-Meier survival curves for each outcome to glyarm
km_fits = lapply(names(outcome_list), function(outcome){
# Extract data of interest
x = outcome_list[[outcome]]
df_sub = na.omit(df[,c(x["time"], x["status"], "glyarm", "kfre_quarts")])
# Formula for survfit object
form = formula(paste0("Surv(", x["time"], ", " ,
x["status"], ") ~ glyarm"))
# Fit survfit object
fit = survfit(form, data = df_sub)
# Manually insert formula
fit$call$formula = form
return(fit)
})
names(km_fits) = names(outcome_list)
# Helper function for formatting labels of the form est (lo, hi)
est_with_ci_text = function(x, est = "Est", lo = "X2.5.", hi = "X97.5.",
digits = 2) {
x = round(x[c(est, lo, hi)], 2)
paste0(x[est], " (", x[lo], ", ", x[hi], ")")
}
invisible(lapply(names(outcome_list), function(outcome){
# Extract data of interest
x = outcome_list[[outcome]]
df_sub = na.omit(df[,c(x["time"], x["status"], "glyarm", "kfre_quarts")])
# Plot cumulative incidence curve with risk table
if (outcome == "egfr40_2consec") {
outcome_title = "ESRD or sustained (consecutive) \neGFR reduction of 40%"
} else {
outcome_title = dat_dict$long[which(dat_dict$short == outcome)]
}
p1 = ggsurvplot(km_fits[[outcome]], data = df_sub, censor = FALSE,
title = outcome_title,
conf.int = TRUE, risk.table = TRUE, fun = "event",
legend.labs = c("Control", "Treatment"))
# Create annotation for hazard ratio and 7-year RMST
p1$plot = p1$plot +
annotate("label", x = 0, y = Inf, vjust = 1, hjust = 0,
size = 4, fontface = 2,
label = paste0("HR = ",
est_with_ci_text(coxph_coeff_tabs[[outcome]][1,-1]),
"\n7-year RMST difference\n = ",
est_with_ci_text(kfre_quart_rmst_boot_tabs[[outcome]][1,-1])))
png(paste0("figs_and_tabs/cum_inc_", outcome, ".png"),
res = 100, width = 600, height = 600)
print(p1)
dev.off()
print(p1)
}))
We plot cumulative incidence curves for each outcome, faceted by quartile.
# Fit Kaplan-Meier survival curves within quartile
km_fits_by_quart = lapply(names(outcome_list), function(outcome){
# Extract data of interest
x = outcome_list[[outcome]]
df_sub = na.omit(df[,c(x["time"], x["status"], "glyarm", "kfre_quarts")])
# Formula for survfit object
form = formula(paste0("Surv(", x["time"], ", " ,
x["status"], ") ~ glyarm"))
all_fits = lapply(1:4, function(i){
# Subset data to quartile
df_sub = df_sub[df_sub$kfre_quarts==i, ]
# Fit survfit object
fit = survfit(form, data = df_sub)
# Manually insert formula
fit$call$formula = form
return(fit)
})
return(all_fits)
})
names(km_fits_by_quart) = names(outcome_list)
invisible(lapply(names(outcome_list), function(outcome){
# Extract data of interest
x = outcome_list[[outcome]]
df_sub = na.omit(df[,c(x["time"], x["status"], "glyarm", "kfre_quarts")])
# Formula for survfit object
form = formula(paste0("Surv(", x["time"], ", " ,
x["status"], ") ~ glyarm"))
all_plots = lapply(1:4, function(i){
# Subset data to quartile
df_sub = df_sub[df_sub$kfre_quarts==i, ]
# Plot cumulative incidence curve with risk table
if (outcome == "Neph2") {
outcome_title = "Development of \nmacro-albuminuria"
} else if (outcome == "Neph3") {
outcome_title = "Renal failure or ESRD \n(dialysis) or SCr>3.3"
} else if (outcome == "egfr40_2consec") {
outcome_title = "ESRD or sustained (consecutive) \neGFR reduction of 40%"
} else {
outcome_title = dat_dict$long[which(dat_dict$short == outcome)]
}
p1 = ggsurvplot(km_fits_by_quart[[outcome]][[i]], data = df_sub, censor = FALSE,
title = paste0(outcome_title, ": \nQuartile ", i),
conf.int = TRUE, risk.table = TRUE, fun = "event",
legend.labs = c("Control", "Treatment"))
# Create annotation for hazard ratio and 7-year RMST
p1$plot = p1$plot +
annotate("label", x = 0, y = Inf, vjust = 1, hjust = 0,
size = 4, fontface = 2,
label = paste0("HR = ",
est_with_ci_text(coxph_coeff_tabs[[outcome]][1+i,-1]),
"\n7-year RMST difference\n = ",
est_with_ci_text(kfre_quart_rmst_boot_tabs[[outcome]][1+i,-1])))
return(p1)
})
# Modify y-axis limits so that they are the same for all plots
y_limits = sapply(all_plots, function(x) {
layer_scales(x$plot)$y$range$range
})
for (i in 1:length(all_plots)) {
all_plots[[i]]$plot = all_plots[[i]]$plot +
ylim(min(y_limits[1,]), max(y_limits[2,]))
}
# Collect all quartile plots for outcome
png(paste0("figs_and_tabs/cum_inc_quart_", outcome, ".png"),
res = 100, width = 1200, height = 1200)
arrange_ggsurvplots(all_plots, print = TRUE,
nrow = 2, ncol = 2)
dev.off()
arrange_ggsurvplots(all_plots, print = TRUE,
nrow = 2, ncol = 2)
}))
# Extract absolute risk for each arm, overall and by quartile
abs_risk_tabs = lapply(names(outcome_list), function(outcome){
out = rbind(
data.frame(Quartile = "Overall",
t(1 - summary(km_fits[[outcome]], times = horizon)$surv)
),
data.frame(Quartile = 1:4,
do.call(rbind,
lapply(km_fits_by_quart[[outcome]], function(x){
t(1 - summary(x, times = horizon)$surv)
}))
))
rownames(out) = NULL
colnames(out) = c("Quartile", "Control", "Treatment")
return(out)
})
names(abs_risk_tabs) = names(outcome_list)
# Function to calculate absolute risk reduction
calc_arr = function(km_fit, time, alpha = 0.05) {
km_fit_summary = summary(km_fit, times = horizon)
est = diff(km_fit_summary$surv)
se = sqrt(sum(km_fit_summary$std.err^2))
return(est + c(0, -1, 1) * qnorm(1 - alpha/2) * se)
}
# Calculate absolute risk reduction, overall and by quartile
arr_risk_tabs = lapply(names(outcome_list), function(outcome){
out = rbind(
data.frame(Quartile = "Overall",
t(calc_arr(km_fits[[outcome]], time = horizon))
),
data.frame(Quartile = 1:4,
do.call(rbind,
lapply(km_fits_by_quart[[outcome]], function(x){
t(calc_arr(x, time = horizon))
}))
))
rownames(out) = NULL
colnames(out) = c("Quartile", "Est", "X2.5.", "X97.5.")
return(out)
})
names(arr_risk_tabs) = names(outcome_list)
# Extract RMST for each arm, overall and by quartile
rmst_by_arm_tabs = lapply(names(outcome_list), function(outcome){
out = rbind(
data.frame(Quartile = "Overall",
Control = overall_rmst_fits[[outcome]]$RMST.arm0$rmst[[1]],
Treatment = overall_rmst_fits[[outcome]]$RMST.arm1$rmst[[1]]
),
data.frame(Quartile = 1:4,
do.call(rbind,
lapply(kfre_quart_rmst_fits[[outcome]], function(x){
t(c(Control = x$RMST.arm0$rmst[[1]],
Treatment = x$RMST.arm1$rmst[[1]]))
}))
))
rownames(out) = NULL
return(out)
})
names(rmst_by_arm_tabs) = names(outcome_list)
# Combine HR, absolute risks, ARR, RMSTs, and RMST difference into table
all_est_tab = lapply(names(outcome_list), function(outcome){
out = data.frame(
Quartile = c("Overall", 1:4),
t(sapply(1:5, function(i) {
c(est_with_ci_text(coxph_coeff_tabs[[outcome]][i,-1]),
t(round(abs_risk_tabs[[outcome]][i,-1], 2)),
est_with_ci_text(arr_risk_tabs[[outcome]][i,-1]),
t(round(rmst_by_arm_tabs[[outcome]][i,-1], 2)),
est_with_ci_text(kfre_quart_rmst_boot_tabs[[outcome]][i,-1]))
}))
)
rownames(out) = NULL
colnames(out) = c("Quartile", "HR (95% CI)",
"Event rate (control)", "Event rate (treatment)",
"ARR (95% CI)",
"RMST (control)", "RMST (treatment)",
"RMST difference (95% CI)")
return(out)
})
names(all_est_tab) = names(outcome_list)
# Format table for printing
est_tab = data.frame(
Outcome = rep(names(all_est_tab),
each = nrow(all_est_tab[[1]])),
do.call(rbind, all_est_tab),
row.names = NULL
)
est_tab = merge(dat_dict[,c("short", "long")], est_tab,
by.x = "short", by.y = "Outcome", all.y = TRUE,
sort = FALSE)[,-1]
names(est_tab) = c("Outcome", "Quartile", "HR (95% CI)",
"Event rate (control)", "Event rate (treatment)",
"ARR (95% CI)",
"RMST (control)", "RMST (treatment)",
"RMST difference (95% CI))")
write.csv(est_tab,
file = "figs_and_tabs/est_tab.csv", row.names = FALSE)
est_tab## Outcome Quartile
## 1 Composite kidney outcome 1
## 2 Composite kidney outcome 2
## 3 Composite kidney outcome 3
## 4 Composite kidney outcome Overall
## 5 Composite kidney outcome 4
## 6 Development of macro-albuminuria Overall
## 7 Development of macro-albuminuria 1
## 8 Development of macro-albuminuria 2
## 9 Development of macro-albuminuria 4
## 10 Development of macro-albuminuria 3
## 11 Renal failure or ESRD (dialysis) or SCr>3.3 Overall
## 12 Renal failure or ESRD (dialysis) or SCr>3.3 1
## 13 Renal failure or ESRD (dialysis) or SCr>3.3 2
## 14 Renal failure or ESRD (dialysis) or SCr>3.3 4
## 15 Renal failure or ESRD (dialysis) or SCr>3.3 3
## 16 All cause death Overall
## 17 All cause death 1
## 18 All cause death 2
## 19 All cause death 3
## 20 All cause death 4
## 21 Cardiovascular death Overall
## 22 Cardiovascular death 2
## 23 Cardiovascular death 3
## 24 Cardiovascular death 4
## 25 Cardiovascular death 1
## 26 1st assisted hypoglycemic event Overall
## 27 1st assisted hypoglycemic event 2
## 28 1st assisted hypoglycemic event 3
## 29 1st assisted hypoglycemic event 4
## 30 1st assisted hypoglycemic event 1
## 31 ESRD or sustained (consecutive) eGFR reduction of 40% Overall
## 32 ESRD or sustained (consecutive) eGFR reduction of 40% 1
## 33 ESRD or sustained (consecutive) eGFR reduction of 40% 2
## 34 ESRD or sustained (consecutive) eGFR reduction of 40% 3
## 35 ESRD or sustained (consecutive) eGFR reduction of 40% 4
## 36 <NA> Overall
## 37 <NA> 1
## 38 <NA> 2
## 39 <NA> 3
## 40 <NA> 4
## HR (95% CI) Event rate (control) Event rate (treatment)
## 1 0.93 (0.63, 1.37) 0.09 0.05
## 2 0.78 (0.59, 1.04) 0.16 0.12
## 3 0.73 (0.54, 0.97) 0.15 0.1
## 4 0.75 (0.65, 0.86) 0.16 0.11
## 5 0.66 (0.53, 0.83) 0.26 0.18
## 6 0.72 (0.61, 0.84) 0.12 0.09
## 7 0.73 (0.39, 1.39) 0.05 0.02
## 8 0.77 (0.56, 1.06) 0.12 0.1
## 9 0.62 (0.49, 0.79) 0.23 0.15
## 10 0.77 (0.56, 1.07) 0.12 0.08
## 11 0.93 (0.73, 1.18) 0.05 0.04
## 12 1.09 (0.68, 1.76) 0.03 0.03
## 13 0.89 (0.52, 1.51) 0.04 0.04
## 14 1.03 (0.69, 1.53) 0.07 0.06
## 15 0.65 (0.38, 1.1) 0.04 0.02
## 16 1.2 (1.04, 1.4) 0.11 0.12
## 17 1.37 (0.91, 2.07) 0.04 0.06
## 18 1.15 (0.84, 1.58) 0.09 0.08
## 19 0.95 (0.7, 1.29) 0.1 0.1
## 20 1.3 (1.03, 1.64) 0.16 0.18
## 21 1.3 (1.05, 1.62) 0.05 0.06
## 22 1.48 (0.89, 2.48) 0.03 0.04
## 23 0.85 (0.55, 1.31) 0.06 0.04
## 24 1.4 (1.01, 1.95) 0.09 0.1
## 25 1.75 (0.88, 3.47) 0.01 0.04
## 26 2.91 (2.56, 3.3) 0.09 0.22
## 27 3.9 (2.85, 5.35) 0.08 0.18
## 28 2.99 (2.32, 3.84) 0.09 0.24
## 29 2.41 (1.97, 2.94) 0.16 0.3
## 30 3.05 (2.22, 4.17) 0.05 0.17
## 31 1.02 (0.9, 1.15) 0.25 0.3
## 32 0.92 (0.7, 1.22) 0.18 0.47
## 33 1.03 (0.81, 1.31) 0.3 0.24
## 34 0.92 (0.73, 1.17) 0.27 0.21
## 35 1.16 (0.93, 1.46) 0.19 0.34
## 36 0.91 (0.81, 1.03) 0.18 0.16
## 37 1.08 (0.77, 1.51) 0.11 0.1
## 38 0.94 (0.73, 1.23) 0.15 0.13
## 39 0.8 (0.62, 1.01) 0.18 0.14
## 40 0.89 (0.73, 1.09) 0.26 0.23
## ARR (95% CI) RMST (control) RMST (treatment)
## 1 0.04 (0, 0.09) 2465.89 2476.03
## 2 0.04 (-0.02, 0.1) 2375.33 2415.39
## 3 0.05 (0, 0.1) 2378.07 2428.15
## 4 0.05 (0.02, 0.07) 2365.05 2413.44
## 5 0.08 (0.03, 0.13) 2211.48 2326.26
## 6 0.04 (0.01, 0.06) 2408.31 2450.8
## 7 0.03 (-0.01, 0.08) 2519.65 2530.64
## 8 0.02 (-0.03, 0.08) 2419.04 2450.36
## 9 0.08 (0.03, 0.12) 2243.62 2360.11
## 10 0.03 (-0.01, 0.08) 2418.51 2451.97
## 11 0.01 (-0.01, 0.02) 2500.38 2505.81
## 12 0 (-0.02, 0.02) 2504.05 2499.95
## 13 0 (-0.04, 0.04) 2509.13 2516.66
## 14 0.01 (-0.02, 0.04) 2486.71 2487.7
## 15 0.02 (0, 0.04) 2503.06 2521.27
## 16 -0.01 (-0.03, 0.01) 2453.08 2429.44
## 17 -0.02 (-0.06, 0.01) 2499.84 2478.81
## 18 0.01 (-0.03, 0.04) 2460.21 2440.68
## 19 0 (-0.04, 0.04) 2443.66 2454.72
## 20 -0.02 (-0.07, 0.03) 2415.23 2358.53
## 21 -0.01 (-0.03, 0.01) 2505.78 2490.36
## 22 -0.01 (-0.02, 0.01) 2519.32 2502.86
## 23 0.01 (-0.02, 0.04) 2493.1 2504.98
## 24 -0.02 (-0.05, 0.02) 2478.8 2443.79
## 25 -0.03 (-0.06, 0.01) 2536.11 2518.29
## 26 -0.13 (-0.15, -0.11) 2427.33 2194.36
## 27 -0.1 (-0.15, -0.05) 2475.67 2259.61
## 28 -0.15 (-0.19, -0.11) 2419.58 2172.5
## 29 -0.15 (-0.19, -0.11) 2334.45 2046.06
## 30 -0.11 (-0.16, -0.07) 2477.46 2305.87
## 31 -0.05 (-0.2, 0.1) 2337.86 2336
## 32 -0.28 (-0.73, 0.16) 2383.63 2403.6
## 33 0.06 (-0.12, 0.23) 2360.03 2355.5
## 34 0.06 (-0.13, 0.24) 2348.4 2361.24
## 35 -0.16 (-0.4, 0.09) 2352.3 2322.99
## 36 0.02 (-0.01, 0.06) 2363.29 2376.76
## 37 0.01 (-0.07, 0.08) 2455.98 2446.84
## 38 0.02 (-0.04, 0.08) 2384.37 2390.91
## 39 0.04 (-0.01, 0.09) 2346.18 2387.56
## 40 0.03 (-0.02, 0.09) 2270.13 2293.4
## RMST difference (95% CI))
## 1 10.14 (-23.9, 43.13)
## 2 40.07 (-4.65, 84.7)
## 3 50.08 (7.95, 92.69)
## 4 48.4 (25.34, 69.6)
## 5 114.77 (58.14, 176.41)
## 6 42.49 (22.52, 61.31)
## 7 10.99 (-8.12, 30.61)
## 8 31.32 (-4.86, 68.63)
## 9 116.49 (62.97, 173.51)
## 10 33.46 (-6.22, 71.98)
## 11 5.43 (-7.02, 16.95)
## 12 -4.1 (-30.36, 19.86)
## 13 7.53 (-14.43, 31.42)
## 14 0.99 (-27.6, 27.66)
## 15 18.22 (-3.14, 39.3)
## 16 -23.64 (-42.2, -6.64)
## 17 -21.03 (-47.7, 6.66)
## 18 -19.53 (-53.51, 12.56)
## 19 11.06 (-20.23, 42.14)
## 20 -56.7 (-100.19, -17.5)
## 21 -15.43 (-28.26, -3.57)
## 22 -16.46 (-40.91, 4.34)
## 23 11.88 (-12.68, 39.15)
## 24 -35.01 (-65.91, -5.41)
## 25 -17.82 (-37.08, 0.4)
## 26 -232.97 (-259.47, -208.03)
## 27 -216.06 (-263.96, -168.89)
## 28 -247.08 (-302.38, -192.26)
## 29 -288.38 (-346.09, -227.82)
## 30 -171.59 (-217.89, -127.53)
## 31 -1.86 (-25.85, 19.87)
## 32 19.97 (-25.65, 66.44)
## 33 -4.53 (-52.13, 41.24)
## 34 12.84 (-30.46, 56.44)
## 35 -29.32 (-75.96, 16.64)
## 36 13.46 (-10.1, 34.64)
## 37 -9.14 (-42.1, 26.45)
## 38 6.54 (-37.73, 48.36)
## 39 41.38 (0.72, 86.08)
## 40 23.27 (-33.49, 77.24)Deciles
We repeat the HTE analysis using KFRE deciles to define the subgroups. The tables summarize the point estimates and 95% bootstrap CIs for 7-year RMST differences and normalized RMST differences by KFRE decile.
# kfre decile thresholds
kfre_dec_thresh = quantile(df$kfre5, seq(0.1, 0.9, by = 0.1))
# Create groups based on thresholds
df$kfre_dec = as.numeric(cut(df$kfre5,
breaks = c(-Inf, as.numeric(kfre_dec_thresh), Inf),
labels = c(1:(length(kfre_dec_thresh)+1))))
# RMST difference by decile
kfre_dec_rmst = sapply(outcome_list, function(x){
df_sub = na.omit(df[,c(x["time"], x["status"], "glyarm", 'kfre_dec')])
if (x["status"] == "egfr40_2consec") {
horizon = 2282
} else {
horizon = 365*7
}
sapply(1:10, function(i){
rmst2(df_sub[df_sub$kfre_dec==i,x["time"]],
df_sub[df_sub$kfre_dec==i,x["status"]],
df_sub[df_sub$kfre_dec==i,"glyarm"],
tau = horizon)$unadjusted.result[1,1]
})
})boot_kfre_dec_rmst = foreach::foreach(
b = 1:B,
.packages = c("survRM2")
) %dopar% {
# Set seed
set.seed(b)
# Create bootstrap sample
boot_idx = sample(nrow(df), replace= TRUE)
df_boot = df[boot_idx,]
# Calculate overall RMST difference
overall_rmst = sapply(outcome_list, function(x){
df_sub = na.omit(df_boot[,c(x["time"], x["status"], "glyarm", "kfre_dec")])
rmst2(df_sub[,x["time"]],
df_sub[,x["status"]],
df_sub[,"glyarm"],
tau = horizon)$unadjusted.result[1,1]
})
# Calculate RMST difference for each decile
kfre_dec_rmst = sapply(outcome_list, function(x){
df_sub = na.omit(df_boot[,c(x["time"], x["status"], "glyarm", "kfre_dec")])
if (x["status"] == "egfr40_2consec") {
horizon = 2282
} else {
horizon = 365*7
}
sapply(1:10, function(i){
rmst2(df_sub[df_sub$kfre_dec==i,x["time"]],
df_sub[df_sub$kfre_dec==i,x["status"]],
df_sub[df_sub$kfre_dec==i,"glyarm"],
tau = horizon)$unadjusted.result[1,1]
})
})
return(list(overall_rmst = overall_rmst,
kfre_dec_rmst = kfre_dec_rmst))
}
save(boot_kfre_dec_rmst, file = "boot_kfre_dec_rmst.rData")# Not normalized
kfre_dec_rmst_boot_ci = lapply(1:10, function(i){
dat = sapply(boot_kfre_dec_rmst, function(x){
x$kfre_dec_rmst[i,]
})
apply(dat, 1, quantile, c(0.025, 0.975))
})
# Make table
kfre_dec_rmst_boot_tabs = lapply(names(outcome_list), function(outcome){
out = rbind(
data.frame(Decile = "Overall",
"Est" = overall_rmst[outcome],
t(kfre_overall_rmst_boot_ci[,outcome])),
data.frame(Decile = 1:10,
"Est" = kfre_dec_rmst[,outcome],
t(sapply(kfre_dec_rmst_boot_ci, function(x){x[,outcome]}))))
rownames(out) = NULL
return(out)
})
names(kfre_dec_rmst_boot_tabs) = names(outcome_list)
# Format table for printing
rmst_dec_tab = data.frame(
Outcome = rep(names(kfre_dec_rmst_boot_tabs),
each = nrow(kfre_dec_rmst_boot_tabs[[1]])),
do.call(rbind, kfre_dec_rmst_boot_tabs),
row.names = NULL
)
rmst_dec_tab = merge(dat_dict[,c("short", "long")], rmst_dec_tab,
by.x = "short", by.y = "Outcome", all.y = TRUE,
sort = FALSE)[,-1]
names(rmst_dec_tab) = c("Outcome", "Decile", "Est", "2.5%", "97.5%")
rmst_dec_tab[,c("Est", "2.5%", "97.5%")] =
round(rmst_dec_tab[,c("Est", "2.5%", "97.5%")], 2)
write.csv(rmst_dec_tab,
file = "figs_and_tabs/rmst_dec_tab.csv", row.names = FALSE)
rmst_dec_tab## Outcome Decile Est
## 1 Composite kidney outcome Overall 48.40
## 2 Composite kidney outcome 3 8.71
## 3 Composite kidney outcome 4 20.11
## 4 Composite kidney outcome 1 35.50
## 5 Composite kidney outcome 2 -23.52
## 6 Composite kidney outcome 7 19.82
## 7 Composite kidney outcome 8 76.14
## 8 Composite kidney outcome 5 74.50
## 9 Composite kidney outcome 6 40.66
## 10 Composite kidney outcome 9 117.08
## 11 Composite kidney outcome 10 166.55
## 12 Development of macro-albuminuria 1 18.89
## 13 Development of macro-albuminuria 2 0.97
## 14 Development of macro-albuminuria Overall 42.49
## 15 Development of macro-albuminuria 5 66.17
## 16 Development of macro-albuminuria 6 21.69
## 17 Development of macro-albuminuria 3 0.69
## 18 Development of macro-albuminuria 4 8.81
## 19 Development of macro-albuminuria 9 123.55
## 20 Development of macro-albuminuria 10 166.73
## 21 Development of macro-albuminuria 7 -2.21
## 22 Development of macro-albuminuria 8 76.28
## 23 Renal failure or ESRD (dialysis) or SCr>3.3 3 12.64
## 24 Renal failure or ESRD (dialysis) or SCr>3.3 4 13.29
## 25 Renal failure or ESRD (dialysis) or SCr>3.3 Overall 5.43
## 26 Renal failure or ESRD (dialysis) or SCr>3.3 1 12.80
## 27 Renal failure or ESRD (dialysis) or SCr>3.3 2 -32.08
## 28 Renal failure or ESRD (dialysis) or SCr>3.3 7 13.57
## 29 Renal failure or ESRD (dialysis) or SCr>3.3 8 -5.17
## 30 Renal failure or ESRD (dialysis) or SCr>3.3 5 1.91
## 31 Renal failure or ESRD (dialysis) or SCr>3.3 6 24.34
## 32 Renal failure or ESRD (dialysis) or SCr>3.3 9 1.01
## 33 Renal failure or ESRD (dialysis) or SCr>3.3 10 10.81
## 34 All cause death 5 -38.36
## 35 All cause death 6 47.71
## 36 All cause death Overall -23.64
## 37 All cause death 1 23.81
## 38 All cause death 2 -57.84
## 39 All cause death 3 -8.49
## 40 All cause death 4 -20.53
## 41 All cause death 9 -42.16
## 42 All cause death 10 -27.86
## 43 All cause death 7 -32.43
## 44 All cause death 8 -59.67
## 45 Cardiovascular death 7 -2.13
## 46 Cardiovascular death 8 -51.38
## 47 Cardiovascular death Overall -15.43
## 48 Cardiovascular death 1 5.14
## 49 Cardiovascular death 2 -28.29
## 50 Cardiovascular death 3 -12.18
## 51 Cardiovascular death 4 -20.96
## 52 Cardiovascular death 5 -28.43
## 53 Cardiovascular death 6 31.45
## 54 Cardiovascular death 9 -23.39
## 55 Cardiovascular death 10 -15.24
## 56 1st assisted hypoglycemic event 9 -255.22
## 57 1st assisted hypoglycemic event 10 -322.17
## 58 1st assisted hypoglycemic event Overall -232.97
## 59 1st assisted hypoglycemic event 1 -156.51
## 60 1st assisted hypoglycemic event 2 -179.27
## 61 1st assisted hypoglycemic event 3 -146.58
## 62 1st assisted hypoglycemic event 4 -246.15
## 63 1st assisted hypoglycemic event 5 -244.96
## 64 1st assisted hypoglycemic event 6 -222.18
## 65 1st assisted hypoglycemic event 7 -223.98
## 66 1st assisted hypoglycemic event 8 -322.58
## 67 ESRD or sustained (consecutive) eGFR reduction of 40% Overall -1.86
## 68 ESRD or sustained (consecutive) eGFR reduction of 40% 1 20.17
## 69 ESRD or sustained (consecutive) eGFR reduction of 40% 2 -31.85
## 70 ESRD or sustained (consecutive) eGFR reduction of 40% 3 38.99
## 71 ESRD or sustained (consecutive) eGFR reduction of 40% 4 -1.02
## 72 ESRD or sustained (consecutive) eGFR reduction of 40% 5 -21.61
## 73 ESRD or sustained (consecutive) eGFR reduction of 40% 6 -2.60
## 74 ESRD or sustained (consecutive) eGFR reduction of 40% 7 7.88
## 75 ESRD or sustained (consecutive) eGFR reduction of 40% 8 -20.27
## 76 ESRD or sustained (consecutive) eGFR reduction of 40% 9 -6.69
## 77 ESRD or sustained (consecutive) eGFR reduction of 40% 10 -10.62
## 78 <NA> 3 7.20
## 79 <NA> Overall 13.46
## 80 <NA> 1 -21.03
## 81 <NA> 2 34.48
## 82 <NA> 7 5.00
## 83 <NA> 4 -11.74
## 84 <NA> 5 -18.58
## 85 <NA> 6 52.00
## 86 <NA> 8 37.16
## 87 <NA> 9 -14.76
## 88 <NA> 10 78.95
## 2.5% 97.5%
## 1 25.34 69.60
## 2 -55.43 68.12
## 3 -59.20 94.31
## 4 -9.60 78.13
## 5 -71.09 24.33
## 6 -45.20 82.04
## 7 -0.95 150.41
## 8 4.43 139.17
## 9 -20.18 108.79
## 10 39.82 199.48
## 11 51.28 295.15
## 12 -2.94 43.91
## 13 -29.93 29.66
## 14 22.52 61.31
## 15 3.37 126.58
## 16 -36.22 82.96
## 17 -47.23 46.09
## 18 -53.19 68.30
## 19 48.33 201.58
## 20 56.66 287.16
## 21 -63.40 57.33
## 22 11.44 143.69
## 23 -30.94 53.54
## 24 -28.24 58.48
## 25 -7.02 16.95
## 26 -26.46 50.40
## 27 -67.33 3.42
## 28 -20.96 45.46
## 29 -46.12 37.26
## 30 -26.69 30.71
## 31 -8.94 60.53
## 32 -27.47 27.56
## 33 -40.83 64.28
## 34 -92.86 7.20
## 35 -2.11 93.33
## 36 -42.20 -6.64
## 37 -8.16 56.33
## 38 -112.47 -2.86
## 39 -56.91 39.78
## 40 -74.17 33.21
## 41 -97.85 11.41
## 42 -99.91 42.92
## 43 -84.23 17.94
## 44 -116.84 0.35
## 45 -37.81 38.23
## 46 -95.68 -3.82
## 47 -28.26 -3.57
## 48 -12.84 24.99
## 49 -67.66 12.53
## 50 -45.90 20.38
## 51 -58.87 16.21
## 52 -64.31 3.39
## 53 -4.53 67.83
## 54 -67.87 16.48
## 55 -69.36 39.69
## 56 -350.61 -162.96
## 57 -427.96 -213.76
## 58 -259.47 -208.03
## 59 -230.74 -94.45
## 60 -246.64 -110.25
## 61 -218.61 -74.88
## 62 -328.18 -166.06
## 63 -316.52 -170.72
## 64 -302.43 -141.34
## 65 -317.99 -138.66
## 66 -408.37 -237.44
## 67 -25.85 19.87
## 68 -25.69 66.34
## 69 -85.17 20.11
## 70 -23.29 93.70
## 71 -63.47 61.74
## 72 -73.14 30.99
## 73 -53.35 57.10
## 74 -51.01 68.09
## 75 -79.31 36.21
## 76 -60.38 51.59
## 77 -74.03 42.28
## 78 -58.05 67.32
## 79 -10.10 34.64
## 80 -66.53 25.80
## 81 -21.42 96.05
## 82 -69.94 83.43
## 83 -81.93 54.57
## 84 -94.39 47.63
## 85 -10.32 113.49
## 86 -41.26 116.84
## 87 -91.99 55.99
## 88 -10.99 171.19# Normalized (subtract overall RMST difference)
kfre_dec_rmst_norm_boot_ci = lapply(1:10, function(i){
dat = sapply(boot_kfre_dec_rmst, function(x){
x$kfre_dec_rmst[i,] - x$overall_rmst
})
apply(dat, 1, quantile, c(0.025, 0.975))
})
# Make a table
kfre_dec_rmst_boot_norm_tabs = lapply(names(outcome_list), function(outcome){
out = rbind(
data.frame(Decile = 1:10,
"Est" = kfre_dec_rmst[,outcome] - overall_rmst[outcome],
t(sapply(kfre_dec_rmst_norm_boot_ci, function(x){x[,outcome]}))))
rownames(out) = NULL
return(out)
})
names(kfre_dec_rmst_boot_norm_tabs) = names(outcome_list)
# Format table for printing
rmst_norm_dec_tab = data.frame(
Outcome = rep(names(kfre_dec_rmst_boot_norm_tabs),
each = nrow(kfre_dec_rmst_boot_norm_tabs[[1]])),
do.call(rbind, kfre_dec_rmst_boot_norm_tabs),
row.names = NULL
)
rmst_norm_dec_tab = merge(dat_dict[,c("short", "long")], rmst_norm_dec_tab,
by.x = "short", by.y = "Outcome", all.y = TRUE,
sort = FALSE)[,-1]
names(rmst_norm_dec_tab) = c("Outcome", "Decile", "Est", "2.5%", "97.5%")
rmst_norm_dec_tab[,c("Est", "2.5%", "97.5%")] =
round(rmst_norm_dec_tab[,c("Est", "2.5%", "97.5%")], 2)
write.csv(rmst_norm_dec_tab,
file = "figs_and_tabs/rmst_norm_dec_tab.csv", row.names = FALSE)
rmst_norm_dec_tab## Outcome Decile Est 2.5%
## 1 Composite kidney outcome 6 -7.74 -68.29
## 2 Composite kidney outcome 7 -28.58 -94.28
## 3 Composite kidney outcome 8 27.74 -41.41
## 4 Composite kidney outcome 1 -12.90 -61.26
## 5 Composite kidney outcome 2 -71.92 -118.99
## 6 Composite kidney outcome 3 -39.69 -107.56
## 7 Composite kidney outcome 4 -28.29 -102.44
## 8 Composite kidney outcome 5 26.10 -41.20
## 9 Composite kidney outcome 10 118.15 11.23
## 10 Composite kidney outcome 9 68.68 -2.81
## 11 Development of macro-albuminuria 10 124.24 18.35
## 12 Development of macro-albuminuria 2 -41.52 -72.46
## 13 Development of macro-albuminuria 9 81.06 12.61
## 14 Development of macro-albuminuria 1 -23.60 -51.80
## 15 Development of macro-albuminuria 6 -20.80 -77.63
## 16 Development of macro-albuminuria 3 -41.80 -89.53
## 17 Development of macro-albuminuria 4 -33.68 -90.49
## 18 Development of macro-albuminuria 5 23.68 -35.90
## 19 Development of macro-albuminuria 7 -44.70 -102.40
## 20 Development of macro-albuminuria 8 33.80 -31.43
## 21 Renal failure or ESRD (dialysis) or SCr>3.3 1 7.36 -31.30
## 22 Renal failure or ESRD (dialysis) or SCr>3.3 3 7.21 -33.70
## 23 Renal failure or ESRD (dialysis) or SCr>3.3 4 7.86 -31.76
## 24 Renal failure or ESRD (dialysis) or SCr>3.3 5 -3.52 -31.38
## 25 Renal failure or ESRD (dialysis) or SCr>3.3 2 -37.51 -72.18
## 26 Renal failure or ESRD (dialysis) or SCr>3.3 7 8.14 -23.65
## 27 Renal failure or ESRD (dialysis) or SCr>3.3 8 -10.60 -48.28
## 28 Renal failure or ESRD (dialysis) or SCr>3.3 9 -4.42 -32.54
## 29 Renal failure or ESRD (dialysis) or SCr>3.3 6 18.91 -13.70
## 30 Renal failure or ESRD (dialysis) or SCr>3.3 10 5.38 -41.88
## 31 All cause death 2 -34.20 -83.89
## 32 All cause death 3 15.15 -33.06
## 33 All cause death 4 3.11 -47.17
## 34 All cause death 1 47.45 12.51
## 35 All cause death 6 71.35 25.75
## 36 All cause death 7 -8.79 -58.82
## 37 All cause death 8 -36.03 -91.25
## 38 All cause death 5 -14.72 -68.10
## 39 All cause death 10 -4.22 -69.41
## 40 All cause death 9 -18.52 -70.92
## 41 Cardiovascular death 6 46.88 11.00
## 42 Cardiovascular death 7 13.30 -23.09
## 43 Cardiovascular death 5 -13.00 -49.17
## 44 Cardiovascular death 1 20.57 0.95
## 45 Cardiovascular death 2 -12.87 -50.65
## 46 Cardiovascular death 3 3.25 -29.87
## 47 Cardiovascular death 4 -5.54 -41.36
## 48 Cardiovascular death 9 -7.96 -48.73
## 49 Cardiovascular death 10 0.19 -48.58
## 50 Cardiovascular death 8 -35.95 -76.10
## 51 1st assisted hypoglycemic event 10 -89.20 -186.05
## 52 1st assisted hypoglycemic event 8 -89.61 -169.00
## 53 1st assisted hypoglycemic event 9 -22.24 -114.53
## 54 1st assisted hypoglycemic event 1 76.46 4.79
## 55 1st assisted hypoglycemic event 2 53.70 -9.64
## 56 1st assisted hypoglycemic event 3 86.39 14.38
## 57 1st assisted hypoglycemic event 4 -13.17 -87.19
## 58 1st assisted hypoglycemic event 5 -11.98 -81.70
## 59 1st assisted hypoglycemic event 6 10.79 -64.52
## 60 1st assisted hypoglycemic event 7 9.00 -79.46
## 61 ESRD or sustained (consecutive) eGFR reduction of 40% 1 22.03 -25.22
## 62 ESRD or sustained (consecutive) eGFR reduction of 40% 2 -29.99 -78.54
## 63 ESRD or sustained (consecutive) eGFR reduction of 40% 3 40.85 -18.43
## 64 ESRD or sustained (consecutive) eGFR reduction of 40% 4 0.84 -59.93
## 65 ESRD or sustained (consecutive) eGFR reduction of 40% 5 -19.75 -70.09
## 66 ESRD or sustained (consecutive) eGFR reduction of 40% 6 -0.74 -49.54
## 67 ESRD or sustained (consecutive) eGFR reduction of 40% 7 9.74 -47.52
## 68 ESRD or sustained (consecutive) eGFR reduction of 40% 8 -18.41 -71.85
## 69 ESRD or sustained (consecutive) eGFR reduction of 40% 9 -4.83 -58.53
## 70 ESRD or sustained (consecutive) eGFR reduction of 40% 10 -8.76 -68.69
## 71 <NA> 6 38.53 -22.03
## 72 <NA> 7 -8.46 -82.39
## 73 <NA> 8 23.69 -46.55
## 74 <NA> 9 -28.23 -99.27
## 75 <NA> 10 65.49 -18.51
## 76 <NA> 2 21.02 -32.91
## 77 <NA> 3 -6.26 -66.91
## 78 <NA> 4 -25.21 -92.36
## 79 <NA> 1 -34.49 -81.03
## 80 <NA> 5 -32.04 -101.18
## 97.5%
## 1 56.03
## 2 32.81
## 3 98.84
## 4 32.11
## 5 -22.45
## 6 15.23
## 7 41.62
## 8 86.52
## 9 232.10
## 10 150.40
## 11 238.20
## 12 -7.63
## 13 155.53
## 14 2.42
## 15 36.38
## 16 4.41
## 17 25.05
## 18 79.58
## 19 12.97
## 20 97.58
## 21 43.90
## 22 46.82
## 23 50.35
## 24 26.20
## 25 -4.79
## 26 37.52
## 27 28.83
## 28 22.22
## 29 53.74
## 30 57.68
## 31 20.44
## 32 62.04
## 33 56.46
## 34 82.73
## 35 115.58
## 36 44.69
## 37 20.33
## 38 27.66
## 39 58.15
## 40 32.59
## 41 84.01
## 42 51.78
## 43 17.93
## 44 42.00
## 45 28.57
## 46 34.99
## 47 33.20
## 48 29.62
## 49 49.38
## 50 10.60
## 51 13.03
## 52 -8.63
## 53 64.96
## 54 139.23
## 55 121.70
## 56 157.80
## 57 63.94
## 58 61.59
## 59 87.80
## 60 99.37
## 61 66.89
## 62 21.12
## 63 96.75
## 64 56.84
## 65 34.63
## 66 57.24
## 67 69.21
## 68 37.15
## 69 49.27
## 70 50.02
## 71 95.67
## 72 66.51
## 73 97.65
## 74 38.38
## 75 149.57
## 76 80.61
## 77 51.53
## 78 37.03
## 79 10.34
## 80 29.27RMST differences.
for (outcome in names(kfre_dec_rmst_boot_tabs)) {
ate_tab = kfre_dec_rmst_boot_tabs[[outcome]][1,]
p1 = kfre_dec_rmst_boot_tabs[[outcome]][-1,] %>%
mutate(x = 1:10) %>%
ggplot(aes(x = x, y = as.numeric(Est))) +
geom_pointrange(aes(ymin = X2.5., ymax = X97.5.), size = 0.25) +
geom_hline(yintercept = ate_tab[,"Est"], color = "grey") +
annotate("rect", xmin = -Inf, xmax = Inf,
ymin = ate_tab[,"X2.5."], ymax = ate_tab[,"X97.5."],
alpha = 0.2) +
scale_x_continuous(breaks = 1:10, labels = 1:10) +
theme(legend.position = 'bottom', legend.box = 'vertical') +
xlab("Decile") + ylab("RMST difference") +
ggtitle(dat_dict$long[which(dat_dict$short == outcome)])
print(p1)
ggsave(paste0("figs_and_tabs/rmst_dec_", outcome, ".png"),
width = 6, height = 4)
}
Normalized RMST differences.
for (outcome in names(kfre_dec_rmst_boot_norm_tabs)) {
p1 = kfre_dec_rmst_boot_norm_tabs[[outcome]] %>%
mutate(x = 1:10) %>%
ggplot(aes(x = x, y = as.numeric(Est))) +
geom_pointrange(aes(ymin = X2.5., ymax = X97.5.), size = 0.25) +
scale_x_continuous(breaks = 1:10, labels = 1:10) +
theme(legend.position = 'bottom', legend.box = 'vertical') +
xlab("Decile") + ylab("Normalized RMST difference") +
geom_hline(yintercept = 0, color = "grey") +
ggtitle(dat_dict$long[which(dat_dict$short == outcome)])
print(p1)
ggsave(paste0("figs_and_tabs/rmst_norm_dec_", outcome, ".png"),
width = 6, height = 4)
}
Jiang’s Diabetic Kidney Disease (DKD) score
We use the score from Jiang (2020) “Establishment and Validation of a Risk Prediction Model for Early Diabetic Kidney Disease Based on a Systematic Review and Meta-Analysis of 20 Cohorts” to define quartile subgroups. There are 1327 missing values, mostly due to missing smoking status.
# Diabetic retinopathy variable
baselinehistoryandphyisical = read.csv("/Users/jliang/Library/CloudStorage/Box-Box/RELATE\ CKD/Study\ Datasets/ACCORD/ACCORD_2017b_2\ 2/Main_Study/4-Data_Sets-CRFs/4a-CRF_Data_Sets/csv/f07_baselinehistoryphysicalexam.csv")
baselinehistoryandphyisical = baselinehistoryandphyisical %>%
mutate(retinopathy = case_when(retpathy==1 | lrtpathy==1 | rrtpathy==1 ~ 1,
TRUE ~ 0))
# Merge in diabetic retinopathy and years of diabetes
df_dkd = merge(df,
baselinehistoryandphyisical %>%
select(MaskID, retinopathy),
by.x = "maskid", by.y = "MaskID", all.x = TRUE)
# Calculate dkd score from Jiang 2020
df_dkd = df_dkd %>%
mutate(dkd = case_when(baseline_age <= 49 ~ 0,
baseline_age <= 59 ~ 3,
baseline_age > 59 ~ 6) +
case_when(bmi < 25 ~ 0,
bmi < 30 ~ 1.5,
bmi >= 30 ~ 3) +
case_when(smokelif == FALSE ~ 0,
smokelif == TRUE ~ 4) +
case_when(retinopathy == FALSE ~ 0,
retinopathy == TRUE ~ 3) +
case_when(hba1c < 7 ~ 0,
hba1c < 8 ~ 1.5,
hba1c < 9 ~ 3,
hba1c >= 9 ~ 4.5) +
case_when(sbp < 130 ~ 0,
sbp < 140 ~ 2,
sbp < 150 ~ 4,
sbp >= 150 ~ 6) +
case_when(hdl < 1.3*18 ~ 0,
hdl >= 1.3*18 ~ 2.5) +
case_when(trig < 1.7*18 ~ 0,
trig >= 1.7*18 ~ 4)+
case_when(uacr < 10 ~ 0,
uacr < 20 ~ 2,
uacr >= 20 ~ 4))
df_dkd %>% summarise_all(~ sum(is.na(.)))## maskid glyarm kfrs neph_composite neph_composite_fu Neph2 Neph2Days Neph3
## 1 0 0 0 1321 1321 1321 1321 0
## Neph3Days ALLDEATH alldeath_fu CVDEATH cvdeath_fu hypoglycemia
## 1 17 0 0 0 0 0
## hypoglycemia_fu cvd_primary cvd_primary_fu female baseline_age raceth sbp dbp
## 1 0 0 0 0 0 0 44 44
## pulsepres chol trig vldl ldl hdl alt cpk fpg hba1c ualb ucreat uacr
## 1 44 2 2 2 2 1 0 1 0 14 0 0 0
## ckd2021GFR screat bmi smokelif anti_htn chol_lowering diuretic insulin
## 1 0 0 6 1266 0 0 0 0
## oral_DM raceth0 raceth1 raceth2 raceth3 days_from_baseline egfr40_2consec
## 1 0 0 0 0 0 0 0
## egfr40_2consec_fu kfre5 kfre_quarts kfre_dec retinopathy dkd
## 1 17 0 0 0 0 1327# Drop NAs
df_dkd = df_dkd %>%
drop_na(dkd)
# DKD quartile thresholds
dkd_quart_thresh = quantile(df_dkd$dkd, seq(0.25, 0.75, by = 0.25))
# Create groups based on thresholds
df_dkd$dkd_quarts = as.numeric(cut(df_dkd$dkd,
breaks = c(-Inf, as.numeric(dkd_quart_thresh), Inf),
labels = c(1:(length(dkd_quart_thresh)+1))))This score seems to be somewhat correlated with 5-year KFRE.
# Spearman's correlation
cor(df_dkd$kfre5, df_dkd$dkd, use = "complete.obs", method = "spearman")## [1] 0.2941996# Scatterplot
df_dkd %>%
ggplot(aes(x = kfre5, y = dkd)) +
geom_point() +
scale_x_log10() +
xlab("5-year predicted risk by KFRE (log)") + ylab("Jiang's DKD score")
ggsave("figs_and_tabs/dkd_vs_kfre.png", width = 6, height = 5)Number of observations overall and in each arm, restricting to those with non-missing covariates needed to calculate the DKD score.
c(table(df_dkd$glyarm), Overall = nrow(df_dkd))## 0 1 Overall
## 4231 4219 8450Number of non-missing observations for each outcome.
t(sapply(outcome_list, function(x) {
my_df = na.omit(df_dkd[,c(x["status"], x["time"], "glyarm")])
c(table(my_df$glyarm), Overall = nrow(my_df))
}))## 0 1 Overall
## neph_composite 3708 3649 7357
## Neph2 3708 3649 7357
## Neph3 4224 4210 8434
## ALLDEATH 4231 4219 8450
## CVDEATH 4231 4219 8450
## hypoglycemia 4231 4219 8450
## cvd_primary 4231 4219 8450
## egfr40_2consec 4224 4210 8434We estimate the 7-year RMST differences for each subgroup defined by the DKD quartiles.
# Overall
dkd_overall_rmst = sapply(outcome_list, function(x){
df_sub = na.omit(df_dkd[,c(x["time"], x["status"], "glyarm", "dkd_quarts")])
rmst2(df_sub[,x["time"]],
df_sub[,x["status"]],
df_sub[,"glyarm"],
tau = horizon)$unadjusted.result[1,1]
})
# By quartile
dkd_quart_rmst = sapply(outcome_list, function(x){
df_sub = na.omit(df_dkd[,c(x["time"], x["status"], "glyarm", "dkd_quarts")])
if (x["status"] == "egfr40_2consec") {
horizon = 2548
} else {
horizon = 365*7
}
sapply(1:4, function(i){
rmst2(df_sub[df_sub$dkd_quarts==i,x["time"]],
df_sub[df_sub$dkd_quarts==i,x["status"]],
df_sub[df_sub$dkd_quarts==i,"glyarm"],
tau = horizon)$unadjusted.result[1,1]
})
})We use 1000 bootstrap samples to obtain confidence intervals for the RMST differences. Tables of the RMST differences and normalized RMST differences (where the overall RMST difference is subtracted from each quartile’s estime) with 95% bootstrap CIs are shown below.
boot_dkd_rmst = foreach::foreach(
b = 1:B,
.packages = c("survRM2")
) %dopar% {
# Set seed
set.seed(b)
# Create bootstrap sample
boot_idx = sample(nrow(df), replace= TRUE)
df_boot = df_dkd[boot_idx,]
# Calculate overall RMST difference
overall_rmst = sapply(outcome_list, function(x){
df_sub = na.omit(df_boot[,c(x["time"], x["status"], "glyarm", "dkd_quarts")])
rmst2(df_sub[,x["time"]],
df_sub[,x["status"]],
df_sub[,"glyarm"],
tau = horizon)$unadjusted.result[1,1]
})
# Calculate RMST difference for each quartile
dkd_quart_rmst = sapply(outcome_list, function(x){
df_sub = na.omit(df_boot[,c(x["time"], x["status"], "glyarm", "dkd_quarts")])
if (x["status"] == "egfr40_2consec") {
horizon = 2548
} else {
horizon = 365*7
}
sapply(1:4, function(i){
rmst2(df_sub[df_sub$dkd_quarts==i,x["time"]],
df_sub[df_sub$dkd_quarts==i,x["status"]],
df_sub[df_sub$dkd_quarts==i,"glyarm"],
tau = horizon)$unadjusted.result[1,1]
})
})
return(list(overall_rmst = overall_rmst,
dkd_quart_rmst = dkd_quart_rmst))
}
save(boot_dkd_rmst, file = "boot_dkd_rmst.rData")# Overall
dkd_overall_rmst_boot_ci = apply(
sapply(boot_dkd_rmst, function(x){
x$overall_rmst
}), 1, quantile, c(0.025, 0.975))
# Not normalized
dkd_quart_rmst_boot_ci = lapply(1:4, function(i){
dat = sapply(boot_dkd_rmst, function(x){
x$dkd_quart_rmst[i,]
})
apply(dat, 1, quantile, c(0.025, 0.975))
})
# Make tables
dkd_quart_rmst_boot_tabs = lapply(names(outcome_list), function(outcome){
out = rbind(
data.frame(Quartile = "Overall",
"Est" = dkd_overall_rmst[outcome],
t(dkd_overall_rmst_boot_ci[,outcome])),
data.frame(Quartile = 1:4,
"Est" = dkd_quart_rmst[,outcome],
t(sapply(dkd_quart_rmst_boot_ci, function(x){x[,outcome]}))))
rownames(out) = NULL
return(out)
})
names(dkd_quart_rmst_boot_tabs) = names(outcome_list)
# Format table for printing
dkd_rmst_tab = data.frame(
Outcome = rep(names(dkd_quart_rmst_boot_tabs),
each = nrow(dkd_quart_rmst_boot_tabs[[1]])),
do.call(rbind, dkd_quart_rmst_boot_tabs),
row.names = NULL
)
dkd_rmst_tab = merge(dat_dict[,c("short", "long")], dkd_rmst_tab,
by.x = "short", by.y = "Outcome", all.y = TRUE,
sort = FALSE)[,-1]
names(dkd_rmst_tab) = c("Outcome", "Quartile", "Est", "2.5%", "97.5%")
dkd_rmst_tab[,c("Est", "2.5%", "97.5%")] =
round(dkd_rmst_tab[,c("Est", "2.5%", "97.5%")], 2)
write.csv(dkd_rmst_tab,
file = "figs_and_tabs/dkd_rmst_tab.csv", row.names = FALSE)
dkd_rmst_tab## Outcome Quartile Est
## 1 Composite kidney outcome 1 -1.45
## 2 Composite kidney outcome 2 85.18
## 3 Composite kidney outcome 3 42.25
## 4 Composite kidney outcome Overall 44.87
## 5 Composite kidney outcome 4 66.30
## 6 Development of macro-albuminuria Overall 39.29
## 7 Development of macro-albuminuria 1 13.15
## 8 Development of macro-albuminuria 2 59.52
## 9 Development of macro-albuminuria 4 74.78
## 10 Development of macro-albuminuria 3 26.68
## 11 Renal failure or ESRD (dialysis) or SCr>3.3 Overall 2.02
## 12 Renal failure or ESRD (dialysis) or SCr>3.3 1 -18.98
## 13 Renal failure or ESRD (dialysis) or SCr>3.3 2 16.80
## 14 Renal failure or ESRD (dialysis) or SCr>3.3 4 -14.72
## 15 Renal failure or ESRD (dialysis) or SCr>3.3 3 23.12
## 16 All cause death Overall -24.05
## 17 All cause death 1 -21.62
## 18 All cause death 2 -13.31
## 19 All cause death 3 -30.15
## 20 All cause death 4 -31.30
## 21 Cardiovascular death Overall -14.34
## 22 Cardiovascular death 2 -13.14
## 23 Cardiovascular death 3 -7.32
## 24 Cardiovascular death 4 -19.76
## 25 Cardiovascular death 1 -17.02
## 26 1st assisted hypoglycemic event Overall -243.40
## 27 1st assisted hypoglycemic event 2 -215.30
## 28 1st assisted hypoglycemic event 3 -265.54
## 29 1st assisted hypoglycemic event 4 -318.35
## 30 1st assisted hypoglycemic event 1 -189.00
## 31 ESRD or sustained (consecutive) eGFR reduction of 40% Overall -6.39
## 32 ESRD or sustained (consecutive) eGFR reduction of 40% 1 -22.44
## 33 ESRD or sustained (consecutive) eGFR reduction of 40% 2 18.58
## 34 ESRD or sustained (consecutive) eGFR reduction of 40% 3 35.63
## 35 ESRD or sustained (consecutive) eGFR reduction of 40% 4 -61.36
## 36 <NA> Overall 14.79
## 37 <NA> 1 2.73
## 38 <NA> 2 21.86
## 39 <NA> 3 -1.66
## 40 <NA> 4 44.36
## 2.5% 97.5%
## 1 -31.25 31.69
## 2 46.13 126.42
## 3 -6.62 98.21
## 4 22.40 70.88
## 5 3.02 135.12
## 6 20.68 61.41
## 7 -9.41 36.59
## 8 26.35 92.45
## 9 14.07 137.76
## 10 -15.88 76.56
## 11 -10.29 15.15
## 12 -41.51 3.25
## 13 -9.23 42.69
## 14 -43.88 16.34
## 15 -3.98 51.16
## 16 -42.16 -5.43
## 17 -46.72 4.35
## 18 -46.95 19.37
## 19 -69.00 7.03
## 20 -78.55 17.71
## 21 -27.16 -0.42
## 22 -41.48 10.65
## 23 -34.33 19.71
## 24 -55.01 17.40
## 25 -35.28 1.49
## 26 -273.12 -215.90
## 27 -264.70 -167.93
## 28 -320.57 -206.13
## 29 -386.64 -249.73
## 30 -237.16 -144.80
## 31 -29.37 18.53
## 32 -62.25 18.13
## 33 -22.70 62.02
## 34 -14.47 87.44
## 35 -122.34 -3.70
## 36 -8.20 38.68
## 37 -33.68 40.86
## 38 -21.75 65.14
## 39 -52.46 46.67
## 40 -18.26 102.03# Normalized (subtract overall RMST difference)
dkd_quart_rmst_norm_boot_ci = lapply(1:4, function(i){
dat = sapply(boot_dkd_rmst, function(x){
x$dkd_quart_rmst[i,] - x$overall_rmst
})
apply(dat, 1, quantile, c(0.025, 0.975))
})
# Make tables
dkd_quart_rmst_boot_norm_tabs = lapply(names(outcome_list), function(outcome){
data.frame(Quartile = 1:4,
"Est" = dkd_quart_rmst[,outcome] - dkd_overall_rmst[outcome],
t(sapply(dkd_quart_rmst_norm_boot_ci, function(x){x[,outcome]})))
})
names(dkd_quart_rmst_boot_norm_tabs) = names(outcome_list)
# Format table for printing
dkd_rmst_norm_tab = data.frame(
Outcome = rep(names(dkd_quart_rmst_boot_norm_tabs),
each = nrow(dkd_quart_rmst_boot_norm_tabs[[1]])),
do.call(rbind, dkd_quart_rmst_boot_norm_tabs),
row.names = NULL
)
dkd_rmst_norm_tab = merge(dat_dict[,c("short", "long")], dkd_rmst_norm_tab,
by.x = "short", by.y = "Outcome", all.y = TRUE,
sort = FALSE)[,-1]
names(dkd_rmst_norm_tab) = c("Outcome", "Quartile", "Est", "2.5%", "97.5%")
dkd_rmst_norm_tab[,c("Est", "2.5%", "97.5%")] =
round(dkd_rmst_norm_tab[,c("Est", "2.5%", "97.5%")], 2)
write.csv(dkd_rmst_norm_tab,
file = "figs_and_tabs/dkd_rmst_norm_tab.csv", row.names = FALSE)
dkd_rmst_norm_tab## Outcome Quartile Est
## 1 Composite kidney outcome 1 -46.32
## 2 Composite kidney outcome 2 40.30
## 3 Composite kidney outcome 3 -2.62
## 4 Composite kidney outcome 4 21.42
## 5 Development of macro-albuminuria 1 -26.14
## 6 Development of macro-albuminuria 2 20.23
## 7 Development of macro-albuminuria 3 -12.61
## 8 Development of macro-albuminuria 4 35.49
## 9 Renal failure or ESRD (dialysis) or SCr>3.3 1 -21.00
## 10 Renal failure or ESRD (dialysis) or SCr>3.3 2 14.78
## 11 Renal failure or ESRD (dialysis) or SCr>3.3 3 21.11
## 12 Renal failure or ESRD (dialysis) or SCr>3.3 4 -16.74
## 13 All cause death 1 2.43
## 14 All cause death 2 10.74
## 15 All cause death 3 -6.11
## 16 All cause death 4 -7.25
## 17 Cardiovascular death 1 -2.68
## 18 Cardiovascular death 2 1.20
## 19 Cardiovascular death 3 7.02
## 20 Cardiovascular death 4 -5.43
## 21 1st assisted hypoglycemic event 1 54.40
## 22 1st assisted hypoglycemic event 2 28.11
## 23 1st assisted hypoglycemic event 3 -22.14
## 24 1st assisted hypoglycemic event 4 -74.95
## 25 ESRD or sustained (consecutive) eGFR reduction of 40% 1 -16.05
## 26 ESRD or sustained (consecutive) eGFR reduction of 40% 2 24.97
## 27 ESRD or sustained (consecutive) eGFR reduction of 40% 3 42.02
## 28 ESRD or sustained (consecutive) eGFR reduction of 40% 4 -54.97
## 29 <NA> 1 -12.06
## 30 <NA> 2 7.07
## 31 <NA> 3 -16.45
## 32 <NA> 4 29.57
## 2.5% 97.5%
## 1 -79.03 -17.43
## 2 3.42 78.62
## 3 -43.95 42.97
## 4 -33.15 77.88
## 5 -53.10 -1.76
## 6 -9.15 50.15
## 7 -49.53 27.74
## 8 -15.36 88.31
## 9 -41.77 -0.67
## 10 -8.62 36.78
## 11 -1.54 44.95
## 12 -41.99 9.48
## 13 -24.89 28.70
## 14 -18.88 37.15
## 15 -38.22 25.78
## 16 -46.67 34.86
## 17 -21.48 17.29
## 18 -22.85 21.44
## 19 -18.10 29.57
## 20 -35.55 27.46
## 21 9.45 94.72
## 22 -14.99 70.82
## 23 -70.53 25.73
## 24 -131.28 -16.13
## 25 -53.80 19.94
## 26 -13.72 59.54
## 27 -1.53 86.23
## 28 -108.37 -5.96
## 29 -46.23 22.40
## 30 -29.95 42.04
## 31 -59.41 23.36
## 32 -21.98 78.94These are plots of the RMST differences by outcome, with a horizontal reference line drawn at the overall/ATE point estimate. Shading denotes the ATE CI.
for (outcome in names(dkd_quart_rmst_boot_tabs)) {
ate_tab = dkd_quart_rmst_boot_tabs[[outcome]][1,]
p1 = dkd_quart_rmst_boot_tabs[[outcome]][-1,] %>%
mutate(x = 1:4) %>%
ggplot(aes(x = x, y = as.numeric(Est))) +
geom_pointrange(aes(ymin = X2.5., ymax = X97.5.), size = 0.25) +
geom_hline(yintercept = ate_tab[,"Est"], color = "grey") +
annotate("rect", xmin = -Inf, xmax = Inf,
ymin = ate_tab[,"X2.5."], ymax = ate_tab[,"X97.5."],
alpha = 0.2) +
scale_x_continuous(breaks = 1:4, labels = 1:4) +
theme(legend.position = 'bottom', legend.box = 'vertical') +
xlab("Quartile") + ylab("RMST difference") +
ggtitle(dat_dict$long[which(dat_dict$short == outcome)])
print(p1)
ggsave(paste0("figs_and_tabs/dkd_rmst_", outcome, ".png"),
width = 6, height = 4)
}
Normalized RMST differences.
for (outcome in names(dkd_quart_rmst_boot_norm_tabs)) {
p1 = dkd_quart_rmst_boot_norm_tabs[[outcome]] %>%
mutate(x = 1:4) %>%
ggplot(aes(x = x, y = as.numeric(Est))) +
geom_pointrange(aes(ymin = X2.5., ymax = X97.5.), size = 0.25) +
scale_x_continuous(breaks = 1:4, labels = 1:4) +
theme(legend.position = 'bottom', legend.box = 'vertical') +
xlab("Quartile") + ylab("Normalized RMST difference") +
geom_hline(yintercept = 0, color = "grey") +
ggtitle(dat_dict$long[which(dat_dict$short == outcome)])
print(p1)
ggsave(paste0("figs_and_tabs/dkd_rmst_norm_", outcome, ".png"),
width = 6, height = 4)
}
